10 Previsiones con modelos de suavizamiento exponencial

En este capítulo, se ampliará el uso de las funciones de suavizado y sus aplicaciones de previsión. La familia de modelos de previsión puede manejar diversos tipos de series temporales, desde las más simples hasta las más complejas. El modelo de medios móviles y los de suavización exponencial se utilizan como punto de partida, y luego se añaden capas al modelo para aumentar su capacidad de manejo de datos complejos de series temporales.

En este capítulo trataremos los siguientes temas: • Predicción con modelos de medias móviles. • Previsiones con modelos de suavizado. • Ajuste de los parámetros de los modelos de suavizado.

#Importar librerias

library(forecast)
## Warning: package 'forecast' was built under R version 4.2.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(h2o)
## Warning: package 'h2o' was built under R version 4.2.3
## 
## ----------------------------------------------------------------------
## 
## Your next step is to start H2O:
##     > h2o.init()
## 
## For H2O package documentation, ask for help:
##     > ??h2o
## 
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit https://docs.h2o.ai
## 
## ----------------------------------------------------------------------
## 
## Attaching package: 'h2o'
## The following objects are masked from 'package:stats':
## 
##     cor, sd, var
## The following objects are masked from 'package:base':
## 
##     %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
##     colnames<-, ifelse, is.character, is.factor, is.numeric, log,
##     log10, log1p, log2, round, signif, trunc
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.2.3
library(plotly)
## Warning: package 'plotly' was built under R version 4.2.3
## 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
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(Quandl)
## Warning: package 'Quandl' was built under R version 4.2.3
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.2.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.2.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last

Previsiones con modelos de medias móviles

Las funciones de media móvil son usadas para suavizar datos de series temporales. Estas pueden utilizarse como modelo de previsión. Se presentarán dos de las funciones de pronóstico de media móvil más comunes: la media móvil simple y la media móvil ponderada.

La media móvil simple

La función de media móvil simple (SMA) para suavizar datos de series temporales puede utilizarse, con unos sencillos pasos, como modelo de previsión. En primer lugar, la estructura de la función SMA para suavizar datos de series temporales es la siguiente:

Grafico

La función de previsión SMA se utiliza cuando la serie de datos no presenta patrones estructurales como componentes de tendencia y estacionales. En este caso, se supone que los valores pronosticados se acercan a las últimas observaciones de la serie.

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
robusta <- Coffee_Prices[,1]
ts_plot(robusta,
        title = "The Robusta Coffee Monthly Prices",
        Ytitle = "Price in USD",
        Xtitle = "Year")

Figura 1. Se observa un ejemplo de la función SMA personalizada.

A continuación, crearemos una función SMA básica utilizando algunas de las funcionalidades del paquete 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)
}

Los argumentos de la función son los siguientes:

• df: La serie de entrada en un formato de marco de datos de dos columnas, donde la primera columna es un objeto Fecha y la segunda son los valores reales de la serie.

• h: El horizonte de la previsión. A efectos del ejemplo siguiente, la función establece las últimas h observaciones como conjunto de prueba. Esto nos permite comparar el rendimiento del modelo.

• m: La longitud de la ventana móvil.

• w: Las ponderaciones de la media, por defecto, utilizando ponderaciones iguales (o media aritmética).

La función sma_forecast tiene los siguientes componentes:

• Tratamiento de errores: Prueba y verifica si los argumentos de entrada de la función son válidos. Si una de las pruebas definidas no es verdadera, detendrá la ejecución de la función y activará un mensaje de error.

• Preparación de datos: Define el objeto data.frame en función de la longitud de la ventana y el horizonte de previsión.

• Cálculo de los datos: Calcula la media móvil simple y devuelve los resultados.

Utilicemos esta función para demostrar el rendimiento de la función SMA. Vamos a pronosticar 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)

Utilizaremos el paquete plotly para trazar los resultados de las diferentes funciones de la media móvil:

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."))

Figura 2 Las principales observaciones del gráfico anterior son las siguientes:

  1. Si la longitud de la ventana móvil es más corta:

• El intervalo de la previsión se aproxima bastante a las observaciones más recientes de la serie.

• Cuanto más rápido converge la previsión a algún valor constante.

2 Si la longitud de la ventana es mayor.

• Más tarda la previsión en converger a un valor constante.

• Puede manejar mejor las perturbaciones y los valores atípicos.

3• Un modelo de previsión SMA con una ventana móvil de longitud 1 es equivalente al modelo de previsión Naïve.

Aunque la función SMA es bastante sencilla de utilizar y de bajo coste de cálculo, tiene algunas limitaciones:

• La capacidad de previsión de la función SMA se limita a un horizonte corto y puede ser deficiente a largo plazo.

• Este método es limitado para datos de series temporales, sin tendencia ni patrones estacionales. Esto afecta principalmente a la media aritmética, que suaviza el patrón estacional y se vuelve plana a largo plazo.

Promedio móvil ponderado

El promedio móvil ponderado (WMA) es una variante del promedio móvil simple (SMA) que utiliza un promedio ponderado en lugar de un promedio aritmético. La principal ventaja del WMA es que permite distribuir el peso de los retrasos en la ventana móvil, lo que puede ser útil cuando la serie de datos presenta alta manipulación con algunos de sus retrasos. La fórmula matemática del WMA se presenta para formalizar su uso.

Grafico

data(USgas)

USgas_df <- ts_to_prophet(USgas)
# El modelo WMA para aplicar todo el peso en el retraso estacional (retraso 1 2):
USgas_fc_m12a <- sma_forecast(USgas_df,
                              h = 24,
                              m = 12,
                              w = c(1, rep(0,11)))
# El modelo WMA para ponderar el primer retraso con O . 2 y el retraso estacional (retraso 12) con 0,8:
USgas_fc_m12b <- sma_forecast(USgas_df,
                              h = 24,
                              m = 12,
                              w = c(0.8, rep(0,10), 0.2))

Se hace uso el paquete plotly para trazar la salida de ambos modelos de WMA:

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"))

Figura 3 En este ejemplo se obtuvo el promedio móvil ponderado, capturando la oscilación estacional de la serie. Se descubrió que al poner toda la ponderación en el desfase estacional se obtiene un modelo Naive,útil para una serie con un patrón estacional dominante.En el segundo ejemplo, ponderamos el promedio entre el retraso más reciente y el retraso estacional.

A pesar de que la media móvil aditiva puede capturar el componente estacional de una serie, no puede capturar la tendencia de la serie debido al efecto de la media. Por lo tanto, este método empezará a perder eficacia una vez que el horizonte de previsión cruce la longitud de la frecuencia de la serie, por ejemplo, más de un año para series mensuales.

Previsión con suavizado exponencial

Los modelos de suavizado exponencial son uno de los enfoques de previsión más populares entre los modelos tradicionales de series temporales. El concepto de suavizado exponencial es similar al de la media móvil, que se basa en la previsión de los valores futuros de la serie promediando las observaciones pasadas. La principal diferencia entre el suavizado exponencial y la media móvil es que el primero promedia todas las observaciones de la serie, mientras que la segunda promedia un subconjunto de m observaciones.

En esta sección nos centraremos en los principales modelos de pronóstico de suavizamiento exponencial:

• Modelo de suavización exponencial simple.

• Modelo de Holt.

• Modelo Holt-Winters.

Modelo de suavizado exponencial simple

El *modelo SES es un modelo de previsión simple de la familia del suavizado exponencial. Este asume que la serie l se mantiene en el mismo nivel y es adecuado para series sin tendencia ni componentes estacionales. Utiliza una media ponderada de las observaciones pasadas para pronosticar los valores futuros de la serie.

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"))

Figura 4. El siguiente ejemplo muestra el decaimiento de los pesos de las observaciones sobre las 15 observaciones más recientes, para unos valores comprendidos entre 0,01 y 1:

Previsión con la función ses

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

• inicial: Define el método para inicializar el valor de Y, que puede calcularse utilizando las primeras observaciones de la serie al establecer el argumento como simple, o estimándolo con el modelo ets (una versión avanzada del modelo Holt-Winters del paquete de previsión) al establecerlo como óptimo.

• alfa: Define el valor del parámetro de suavizado del modelo. Si se establece en NULL, la función lo estimará.

• h: Define el horizonte de previsión.

Utilicemos de nuevo la función ses para pronosticar los precios mensuales del café Robusta. Dejaremos los últimos 12 meses de la serie como conjunto de prueba para evaluar el rendimiento del modelo. Para ello, utilizaremos la función ts_split del paquete TSstudio:

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

library(forecast)
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

Podemos revisar el rendimiento del modelo usando la función test_forecast:

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)))

Figura 5 Como podemos ver en el gráfico de previsión anterior, la función ses está utilizando el conjunto de entrenamiento para identificar el nivel de la serie mediante la estimación del parámetro alfa y el nivel inicial del modelo (o Ý1). El nivel del valor pronosticado se aproxima bastante al valor de la última observación de la serie, ya que el valor α, en este caso, se aproxima a 1. Dado que el objetivo del modelo ses es pronosticar el nivel de la serie, el modelo no captará ninguna oscilación a corto plazo.

En el caso de una previsión plana, los intervalos de confianza del modelo desempeñan un papel fundamental, ya que el nivel de incertidumbre es mayor. Por lo tanto, será útil evaluar si los valores de previsión están dentro de los límites de los intervalos de confianza del modelo. Utilizaremos la función plot_forecast del paquete TSstudio para crear un gráfico interactivo para el modelo fs_ses que hemos creado y trazar el conjunto de pruebas:

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)))

Figura 6. Como puede verse en el gráfico de previsiones anterior, el conjunto de pruebas se encuentra dentro del intervalo de confianza del 80%.

Optimización de modelos con búsqueda en cuadrícula

La función SES optimiza los valores de los parámetros (α y l0) del modelo que minimizan la suma de errores cuadráticos en el conjunto de entrenamiento. En el caso del modelo SES, también se aplicaría una búsqueda en cuadrícula para identificar el valor óptimo de la constante de suavizado que minimice alguna métrica de error del modelo.

En el ejemplo mencionado, se utiliza una búsqueda en cuadrícula para ajustar los parámetros α y l0 del modelo SES para minimizar el error de la serie de precios de Robusta. Se divide el modelo en particiones de entrenamiento, prueba y validación, y se evalúa el rendimiento del modelo para cada valor de α en la partición de prueba en base a los siguientes pasos:

  1. Entrenar el modelo SES utilizando la partición de entrenamiento y pronosticar las observaciones de la partición de prueba.

  2. Evaluar el rendimiento del modelo tanto en la partición de entrenamiento como en la de pruebas.

  3. Añada las particiones de entrenamiento y de prueba (en orden cronológico) y vuelva a entrenar el modelo en la nueva partición (entrenamiento y prueba) antes de pronosticar los valores de la partición de validación.

  4. Evaluar el rendimiento del segundo modelo en la partición de validación.

Estos pasos nos permitirán seleccionar el valor alfa que minimiza el MAPE en la partición de prueba y utilizar el conjunto de validación para validar el rendimiento del modelo. Antes de establecer la función, establezcamos las particiones de entrenamiento, prueba y validación utilizando la función ts_split:

robusta_par1 <- ts_split(robusta, sample.out = 24)

train1 <- robusta_par1$train 
test1 <- ts_split(robusta_par1$test, sample.out = 12)$train

robusta_par2 <- ts_split(robusta, sample.out = 12)

train2 <- robusta_par2$train
valid <- robusta_par2$test

Utilizaremos las variables train1 y test 1 para la partición de entrenamiento y prueba, y train2 para volver a entrenar el modelo y validar los resultados en la partición válida. La siguiente variable alfa define el rango de búsqueda. Asignaremos una secuencia de valores entre 0 y 1 con un incremento de 0. 01 utilizando la función seq:

alpha <- seq(from = 0, to = 1, by = 0.01)

Como el valor de alfa debe ser mayor que cero, sustituiremos 0 por un número pequeño que se acerque bastante a cero:

alpha[1] <- 0.001

Utilizaremos la función lapply para iterar sobre el modelo utilizando los diferentes valores de a:

ses_grid <- lapply(alpha, function(i){
  md1 <- md_accuracy1 <- md2 <- md_accuracy2 <- results <-  NULL
  md1 <- ses(train1, h = 12, alpha = i, initial = "simple")
  md_accuracy1 <- accuracy(md1, test1)
  
  md2 <- ses(train2, h = 12, alpha = i, initial = "simple")
  md_accuracy2 <- accuracy(md2, valid)
  
  resutls <- data.frame(alpha = i, 
                        train = md_accuracy1[9], 
                        test = md_accuracy1[10], 
                        valid = md_accuracy2[10])
  
}) %>% bind_rows()

Como se puede ver en los siguientes resultados de las pruebas, mientras que α = 1 minimiza el MAPE en la partición de entrenamiento, α = 0,03 minimiza la tasa de error en la partición de pruebas partición de entrenamiento, a = 0,03 minimiza la tasa de error en la partición de prueba. Los resultados en la partición de validación siguen el mismo patrón que la partición de prueba, con una puntuación MAPE del 9,98% en la partición de prueba y del 6,60% en la partición de validación:

plot_ly(data = ses_grid, x = ~ alpha, y = ~ train, 
        line = list(color = 'rgb(205, 12, 24)'),
        type = "scatter", 
        mode = "lines", 
        name = "Training") %>%
  add_lines(x = ~ alpha, y = ~ test, line = list(color = "rgb(22, 96, 167)", dash = "dash"), name=  "Testing") %>%
  add_lines(x = ~ alpha, y = ~ valid, line = list(color = "green", dash = "dot"), name = "Validation") %>%
  layout(title = "SES Model Grid Search Results",
         yaxis = list(title = "MAPE (%)"))

Figura 7

Método Holt

El método Holt, también conocido como modelo de suavización exponencial doble, es una versión ampliada del modelo SES. Se basa en la estimación del nivel y la tendencia más recientes mediante dos parámetros de suavización, α y β. Una vez que el modelo estima el nivel y la tendencia más recientes (LT, y TT, respectivamente), los utiliza para construir la previsión de la serie mediante la siguiente ecuación:

• Para una serie con tendencia aditiva: Grafico

• Para una serie con tendencia multiplicativa: Grafico

Previsión con la función holt

El paquete de forecast proporciona una implementación del modelo Holt con la función holt. Esta función inicializa automáticamente los valores de ´L1 y ´T1, e identifica los valores de α y β que minimizan el SSE del modelo en el conjunto de entrenamiento. En el siguiente ejemplo, recuperaremos los datos trimestrales del Producto Interior Bruto (PIB) de EE.UU. de la API de Datos Económicos de la Reserva Federal (FRED) utilizando el paquete **Quandl*:

library(Quandl)

gdp <- Quandl("FRED/GDP", start_date = "2010-01-01", type = "ts")

ts_info(gdp)
##  The gdp series is a ts object with 1 variable and 48 observations
##  Frequency: 4 
##  Start time: 2010 1 
##  End time: 2021 4
ts_plot(gdp, 
        title = "Gross Domestic Product",
        Ytitle = "Billions of Dollars",
        Xtitle = "Source: U.S. Bureau of Economic Analysis / fred.stlouisfed.org")

Figura 8. Observará en el siguiente gráfico que la serie del PIB tiene una fuerte tendencia lineal y ningún componente estacional (ya que la serie está desestacionalizada):

Como hicimos anteriormente, dejaremos los últimos ocho trimestres para probar y entrenar el modelo con el resto de las observaciones de la serie con la función ho lt:

gdp_par <- ts_split(gdp, sample.out = 8)

train <- gdp_par$train
test <- gdp_par$test
fc_holt <- holt(train,  h = 8, initial = "optimal") 

fc_holt$model
## Holt's method 
## 
## Call:
##  holt(y = train, h = 8, initial = "optimal") 
## 
##   Smoothing parameters:
##     alpha = 0.9985 
##     beta  = 0.0941 
## 
##   Initial states:
##     l = 14603.0487 
##     b = 161.9411 
## 
##   sigma:  80.9869
## 
##      AIC     AICc      BIC 
## 504.8838 506.6485 513.3282

Comparemos el rendimiento del modelo en las particiones de entrenamiento y prueba con la función accuracy:

accuracy(fc_holt, test)
##                      ME       RMSE      MAE         MPE      MAPE       MASE
## Training set   10.99164   76.83093  61.7626  0.05241543 0.3460477 0.08790283
## Test set     -665.87355 1141.19097 854.8069 -3.29699460 4.0874667 1.21659297
##                   ACF1 Theil's U
## Training set 0.0192117        NA
## Test set     0.2975704   1.08055

Los cambios en el crecimiento de la tendencia y la previsión pueden observarse con la función test_forecast:

test_forecast(gdp, forecast.obj = fc_holt, test = test)

Figura 9 Aunque el modelo Holt se diseñó para manejar series temporales con tendencia lineal, el argumento exponencial de la función holt ofrece la opción de manejar series con tendencias exponenciales o decrecientes cuando se establece en TRUE. En el ejemplo anterior, podemos utilizar el argumento exponencial para modificar el patrón de crecimiento de la tendencia.

En este caso, nos gustaría tener un peso más alto para la tendencia, y fijaremos β en 0,75 (un enfoque más sólido para identificar el valor óptimo de β sería utilizar una búsqueda de cuadrícula):

fc_holt_exp <- holt(train,
                    h = 8,
                    beta = 0.75 ,
                    initial = "optimal",
                    exponential = TRUE) 

fc_holt_exp$model
## Holt's method with exponential trend 
## 
## Call:
##  holt(y = train, h = 8, initial = "optimal", exponential = TRUE,  
## 
##  Call:
##      beta = 0.75) 
## 
##   Smoothing parameters:
##     alpha = 0.781 
##     beta  = 0.75 
## 
##   Initial states:
##     l = 14572.2291 
##     b = 1.0142 
## 
##   sigma:  0.0053
## 
##      AIC     AICc      BIC 
## 515.9291 517.0720 522.6847
accuracy(fc_holt_exp, test)
##                       ME       RMSE       MAE         MPE      MAPE      MASE
## Training set   -2.792636   90.20753  74.14323 -0.01756593 0.4165413 0.1055234
## Test set     -703.014806 1151.21262 858.97402 -3.46212317 4.1130726 1.2225238
##                    ACF1 Theil's U
## Training set -0.1541725        NA
## Test set      0.2834645  1.090637

Similarmente, podemos trazar los valores ajustados y previstos frente a los datos reales con las funciones de test_forecast.

test_forecast(gdp, forecast.obj = fc_holt, test = test)

Figura 10 Como puede verse, el índice de error del segundo modelo Holt es más equilibrado, ya que la relación entre el error en el conjunto de pruebas y el de entrenamiento es de 2,5 y 2,1 para las métricas RMSE y MAPE, respectivamente.

El modelo de Holt-Winters

Cerraremos este capítulo con el tercer modelo, y el más avanzado, de la familia de modelos de previsión de suavizado exponencial: el modelo de Holt-Winters. El smodelo de Holt-Winter (HW) es una versión ampliada del modelo de Holt y puede manejar datos de series temporales con componentes de tendencia y estacionales. La previsión del componente estacional requiere un tercer parámetro y ecuación de suavizado, además de los del nivel y la tendencia.

Ambos componentes, tendencial y estacional, pueden tener una estructura aditiva o multiplicidad, lo que añade cierta complejidad al modelo, ya que existen múltiples combinaciones posibles. combinaciones:

• Tendencia aditiva y componentes estacionales.

• Tendencia aditiva y componentes estacionales multiplicativos.

• Tendencia multiplicativa y componentes estacionales aditivos.

• Tendencia multiplicativa y componentes estacionales.

Usemos la función de descomposición para diagnosticar la estructura de la tendencia y los componentes estacionales de la serie:

data(USgas)

decompose(USgas) %>% plot()

Figura 11 Podemos observar en el gráfico anterior que tanto la tendencia como los componentes estacionales de la serie tienen una estructura aditiva.

Al igual que hicimos anteriormente, crearemos particiones de entrenamiento y prueba utilizando los últimos 12 meses de la serie para evaluar el rendimiento del modelo:

USgas_par <- ts_split(USgas, 12)

train <- USgas_par$train
test <- USgas_par$test

A continuación, utilizaremos el modelo Holtwinters para pronosticar los últimos 12 meses de la serie para el conjunto de pruebas):

md_hw <- HoltWinters(train)

md_hw
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = train)
## 
## Smoothing parameters:
##  alpha: 0.371213
##  beta : 0
##  gamma: 0.4422456
## 
## Coefficients:
##             [,1]
## a   2491.9930104
## b     -0.1287005
## s1    32.0972651
## s2   597.1088003
## s3   834.9628439
## s4   353.8593860
## s5   318.1927058
## s6  -173.0721496
## s7  -346.6229990
## s8  -329.7169608
## s9  -112.1664217
## s10 -140.3186476
## s11 -324.5343787
## s12 -243.9334551

El siguiente paso consiste en pronosticar los próximos 12 meses o los valores del conjunto de prueba y evaluar el rendimiento de los modelos con las funciones accuracy y forecast.

fc_hw <- forecast(md_hw, h = 12)

accuracy(fc_hw, test)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  7.828361 115.2062 87.67420 0.2991952 4.248131 0.7614222
## Test set     51.013877 115.1555 98.06531 1.7994297 3.766099 0.8516656
##                     ACF1 Theil's U
## Training set  0.21911103        NA
## Test set     -0.01991923 0.3652142
test_forecast(actual = USgas, 
              forecast.obj = fc_hw, 
              test = test)

Figura 12 El trazado de los valores ajustados y previstos nos permite comprender mejor el rendimiento del modelo. Como podemos ver en el gráfico anterior, el modelo HW captura bien los patrones estacionales de ambas series. Por otra parte, el modelo no capta el pico del año durante el mes de enero en la mayoría de los casos.

Alternativamente, el modelo se puede entrenar con una búsqueda de cuadrícula de forma similar a lo que hicimos con el modelo SES. En este caso, hay tres parámetros que optimizar: α, β y γ. El paquete TSstudio proporciona una función de búsqueda en cuadrícula personalizada basada en el enfoque de backtesting para entrenar una función HoltWinters.

Por razones de eficacia, empezaremos con una búsqueda poco profunda que incluya incrementos mayores en los valores de los parámetros. Esto nos ayudará a acotar el área de búsqueda y luego aplicar una búsqueda más profunda en esas áreas. La búsqueda superficial incluirá backtesting sobre 6 periodos diferentes utilizando una secuencia entre o y 1 con un incremento de 0. 1:

shallow_grid <-  ts_grid(train, 
                         model = "HoltWinters",
                         periods = 6,
                         window_space = 6,
                         window_test = 12,
                         hyper_params = list(alpha = seq(0,1,0.1),
                                             beta = seq(0,1,0.1),
                                             gamma = seq(0,1,0.1)),
                         parallel = TRUE,
                         n.cores = 8)
## Warning in ts_grid(train, model = "HoltWinters", periods = 6, window_space = 6,
## : The value of the 'alpha' parameter cannot be equal to 0 replacing 0 with 1e-5
## Warning in ts_grid(train, model = "HoltWinters", periods = 6, window_space = 6,
## : The value of the 'beta' parameter cannot be equal to 0 replacing 0 with 1e-5
## Warning in ts_grid(train, model = "HoltWinters", periods = 6, window_space = 6,
## : The value of the 'gamma' parameter cannot be equal to 0 replacing 0 with 1e-5
## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.
## Warning: Detected creation of a 'multiprocess' future, which are defunct in
## future (>= 1.32.0) [2023-03-06]. Instead, specify either 'multisession'
## (recommended) or 'multicore'. If still used, 'multiprocess' becomes the same as
## 'sequential'.
## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.

La salida de la rejilla siguiente proporciona cualquier combinación de α, β y Y con la tasa de error en cada período de prueba y su media global. La tabla ordena la media global del modelo de la mejor combinación a la peor:

shallow_grid$grid_df[1:10,] 
##    alpha  beta gamma        1        2        3        4        5        6
## 1    0.7 1e-05   0.9 2.621397 2.184892 4.338268 4.401508 8.634681 5.977068
## 2    0.7 1e-05   1.0 3.132417 2.394857 4.486892 4.080787 7.887104 6.536839
## 3    0.7 1e-05   0.8 2.582036 2.352049 4.543924 4.709976 9.453191 5.425970
## 4    0.7 1e-01   0.9 3.218986 2.808945 4.906638 4.149888 9.431218 4.723020
## 5    0.6 1e-01   0.6 4.077964 1.936141 4.374585 4.114277 9.710363 5.054511
## 6    0.6 1e-05   0.6 2.766628 2.636940 4.747228 4.420048 8.803412 6.386971
## 7    0.6 1e-05   0.7 3.612994 2.990184 4.723011 4.084390 8.023247 6.936358
## 8    0.6 1e-05   0.5 2.581408 2.653861 5.054216 4.765605 9.763453 5.813665
## 9    0.7 1e-01   1.0 4.439969 1.933397 6.896359 4.247895 7.638064 5.918909
## 10   0.6 2e-01   0.6 4.399195 3.011563 5.440750 3.986332 9.271275 5.039445
##        mean
## 1  4.692969
## 2  4.753149
## 3  4.844524
## 4  4.873116
## 5  4.877974
## 6  4.960205
## 7  5.061697
## 8  5.105368
## 9  5.179099
## 10 5.191427

El plot_grid proporciona una visión intuitiva del rango óptimo de valores de cada parámetro mediante un gráfico de paracordes. Por defecto, la función resalta el 10% superior modelos:

plot_grid(shallow_grid)
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 11 > 1' in coercion to
## 'logical(1)'

## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 11 > 1' in coercion to
## 'logical(1)'

## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 11 > 1' in coercion to
## 'logical(1)'

Figura 13 Observará en el gráfico de búsqueda que el intervalo óptimo de a varía entre 0,1 y 0,5, β entre 0 y 0,1, y Y entre 0,2 y 0,4. Esto nos ayudará a establecer el intervalo de búsqueda de hiperparámetros para una búsqueda más profunda. Esto nos ayudará a establecer el rango de búsqueda de hiperparámetros para la búsqueda de cuadrícula más profunda estrechando el rango de búsqueda pero utilizando una búsqueda más granular:

deep_grid <-  ts_grid(train, 
                      model = "HoltWinters",
                      periods = 6,
                      window_space = 6,
                      window_test = 12,
                      hyper_params = list(alpha = seq(0.1,0.5,0.01),
                                          beta = seq(0,0.1,0.01),
                                          gamma = seq(0.2,0.4,0.01)),
                      parallel = TRUE,
                      n.cores = 8)
## Warning in ts_grid(train, model = "HoltWinters", periods = 6, window_space = 6,
## : The value of the 'beta' parameter cannot be equal to 0 replacing 0 with 1e-5
## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.
## Warning: Detected creation of a 'multiprocess' future, which are defunct in
## future (>= 1.32.0) [2023-03-06]. Instead, specify either 'multisession'
## (recommended) or 'multicore'. If still used, 'multiprocess' becomes the same as
## 'sequential'.
## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.
plot_grid(deep_grid)
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 41 > 1' in coercion to
## 'logical(1)'
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 11 > 1' in coercion to
## 'logical(1)'
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 21 > 1' in coercion to
## 'logical(1)'

Figura 14 Como puede verse en este gráfico, el rango de error de los modelos del 10% superior ha disminuido con respecto al de la búsqueda superficial:

plot_grid(deep_grid, type = "3D", top = 250)
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 41 > 1' in coercion to
## 'logical(1)'
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 11 > 1' in coercion to
## 'logical(1)'
## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 21 > 1' in coercion to
## 'logical(1)'

Figura 15 Este gráfico puede visualizarse en una vista 3D (cuando se realiza una búsqueda de tres parámetros):

El último paso de este proceso consiste en extraer los valores de los parámetros de suavizado óptimos del modelo de cuadrícula basado en la búsqueda, volver a entrenar el modelo HW y utilizarlo para pronosticar los valores futuros de las series.

md_hw_grid <- HoltWinters(train, 
                          alpha = deep_grid$alpha,
                          beta = deep_grid$beta,
                          gamma = deep_grid$gamma)

fc_hw_grid <- forecast(md_hw_grid, h = 12)

accuracy(fc_hw_grid, test)
##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 4.106435 117.0347 90.37810  0.2004600 4.387451 0.7849046 0.1526260
## Test set     3.124700 113.5063 97.20942 -0.2685092 3.899686 0.8442324 0.1506711
##              Theil's U
## Training set        NA
## Test set     0.3687748
test_forecast(actual = USgas, 
              forecast.obj = fc_hw_grid, 
              test = test)

**Figura 17* Como puede observarse en los resultados de las funciones accuracy y test_forecast, la optimización del modelo mediante la búsqueda en cuadrícula proporciona, en este caso, una mejora del rendimiento del modelo al reducir la puntuación MAPE del 6,95% al 5,92%.

Resumen

En este capítulo, hemos introducido el uso de una media ponderada de observaciones pasadas para pronosticar datos de series temporales. Comenzamos con un enfoque de previsión simplista e ingenuo con la función de media móvil. Aunque esta función se limita a previsiones a corto plazo y sólo puede manejar series temporales sin componentes estacionales ni de tendencia, proporciona contexto para las funciones de suavizado exponencial. La familia de modelos de previsión del alisamiento exponencial se basa en el uso de distintos parámetros de alisamiento, es decir, o, , y , para modelar los principales componentes de las series temporales: nivel de datos, tendencia y estacional, respectivamente. Las principales ventajas de las funciones de suavizado exponencial son su sencillez, su bajo coste computacional y su modularidad, que les permite manejar distintos tipos de datos de series temporales, como tendencias lineales y exponenciales y componentes estacionales.