PROY 3 - 10

Pronóstico con Exponencial

Modelos de suavizado

En este capitulo se ampliara el uso de las funciones de suavizado y se presentara sus aplicaciones de pronostico. Los modelos de pronostico manejan una gran variedad de tipos de series de tiempo, por lo que comenzaremos con el modelo de promedio movil basico y los modelos de siavizado exponencial simple y luego se va ampliando la dificultad.

En este capitulo se trataran estos temas: -Pronóstico con modelos de promedio móvil -Enfoques de pronóstico con modelos suavizantes -Parámetros de ajuste para modelos de suavizado

Importar librerias

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

Pronóstico con modelos de promedio móvil

La media móvil simple

La función de promedio móvil simple (SMA), que usamos en el Capítulo 5, Descomposición de datos de series de tiempo, para obtener datos de series de tiempo uniformes se puede utilizar, con algunos pasos simples, como modelo de pronóstico.

Se recomienda el pronóstico con la función SMA cuando la serie de entrada no tiene patrones estructurales, como componentes de tendencia y estacionales. En este caso, es razonable suponer que los valores pronosticados están relativamente cerca de 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. En esta se usa una función SMA personalizada y para pronosticar los precios mensuales del café Robusta en los próximos 12 meses. Los precios del café Robusta (USD por kg) son un ejemplo de datos de series de tiempo que no tienen una tendencia específica o patrones estacionales. Más bien, esta serie tiene un componente de ciclo, donde la magnitud y la duración del ciclo siguen cambiando de un ciclo a otro. Esta serie es parte del conjunto de datos Coffee_Pri ces y está disponible en el paquete TSstudio

Ahora se creara una función SMA básica utilizando algunas de las funciones 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 del pronóstico. A los efectos del siguiente ejemplo, la función establece las últimas h observaciones como un conjunto de prueba. Esto nos permite comparar el rendimiento del modelo.

-m: La longitud de la ventana rodante.

-w: Los pesos del promedio, por defecto, usando pesos iguales (o promedio aritmético).

La función sma_forecast tiene los siguientes componentes:

0 Manejo de errores: Pruebe y verifique 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 generará un mensaje de error.

0 Preparación de datos: Esto define los datos. objeto de marco basado en la longitud de la ventana y el horizonte de pronóstico.

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

Utilicemos esta función para demostrar el rendimiento de la función SMA. Pronosticaremos los últimos 2 4 meses de la serie Robusta utilizando una ventana móvil de 3, 6, 12, 24 y 3 6 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)

Hacemos uso de la libreria Plotly:

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. De la anterior figura se observa: Si la longitud de la ventana móvil es más corta: -El rango del pronóstico es bastante cercano a las observaciones más recientes de la serie -Cuanto más rápido converge el pronóstico a algún valor constante c Si la longitud de la ventana es mayor: -Cuanto más tarde hasta que el pronóstico converja a algún valor constante -Puede manejar mejor choques y 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.

Si bien la función SMA es bastante simple de usar y económica en potencia de cómputo, tiene algunas limitaciones: -El poder de pronóstico de la función SMA se limita a un horizonte corto y puede tener un rendimiento deficiente a largo plazo. -Este método está limitado para datos de series de tiempo, sin tendencias ni patrones estacionales. Esto afecta principalmente al promedio aritmético que suaviza el patrón estacional y se vuelve plano a largo plazo.

Media móvil ponderada

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 oposición al 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 rezagos. En este caso, utilizaremos el argumento in para establecer el peso promedio y, por lo tanto, transformar la función de SMA a WMA. Como hicimos w anteriormente, primero transformaremos la serie en datos. formato de cuadro con la función ts_to_prophet:

data(USgas)

USgas_df <- ts_to_prophet(USgas)

Ahora, podremos elegir alguna de las dos estrategias:

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

Otra vez hacemos uso de la libreria Plotly (para trazar la salida de ambos modelos 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. Como se observa, ambos modelos capturaron en cierta medida la oscilación estacional de la serie. Establecer el peso completo 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 el gas estadounidense. En el segundo ejemplo, ponderamos el promedio entre el retraso más reciente y el retraso estacional. Si bien 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).

Pronóstico con suavizado exponencial

En los modelos de pronóstico de series de tiempo tradicionales, la función de suavizado exponencial es uno de los métodos de pronóstico más populares. Este método es conceptualmente similar al método de promedio móvil que presentamos anteriormente, ya que ambos se basan en pronosticar valores futuros de una serie promediando observaciones pasadas de la serie. La principal diferencia entre el suavizado exponencial y el promedio móvil es que el primero toma el promedio de todas las observaciones de la serie, en lugar del segundo un subconjunto de m observaciones.

En esta sección nos centraremos en los principales modelos de pronóstico de suavizamiento exponencial: -Modelo de suavizado exponencial simple -Modelo Holt -Modelo Holt Inviernos.

Modelo de suavizado exponencial simple

El suavizado exponencial simple (SES), como su nombre lo indica, es el modelo de pronóstico más simple entre la familia de suavizado exponencial. El supuesto principal 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.

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 demuestra el decaimiento de los pesos de las observaciones en las 1 5 observaciones más recientes, para 04 valores entre 0 . 0 1 a 1.

Predicción con la función ses

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: -Inicial: define el método para inicializar el valor de Y1, que se puede calcular utilizando las primeras observaciones de la serie estableciendo el argumento en simple, o estimándolo con el modelo ets (una versión avanzada de Holt-Winters). modelo del paquete de pronóstico) al establecerlo en óptimo. -Alfa: 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 previsión.

Usemos 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 comparar el rendimiento del modelo. Haremos esto usando 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 pronóstico anterior, la función ses utiliza el conjunto de entrenamiento para identificar el nivel de la serie al estimar el parámetro alfa y el nivel inicial del modelo (o Y1). El nivel de valor de pronóstico es bastante cercano al valor de la última observación de la serie ya que el valor (1’, en este caso, es cercano a 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.

En el caso de un pronóstico plano, el intervalo de confianza del modelo juega un papel crucial porque el nivel de incertidumbre es mayor. Por lo tanto, es útil evaluar si los valores pronosticados están dentro del intervalo de confianza del modelo. Usaremos la función plot_forecast del paquete TSstudio para crear gráficos interactivos y trazar el conjunto de pruebas para el modelo f s_ses que creamos:

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. Si observamos, notaremos que el conjunto de prueba está dentro del rango 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 del modelo (oz y lo) que minimizan la suma de los errores cuadráticos (SSE) del modelo en el conjunto de entrenamiento. En el caso del modelo SES, aplicaremos una búsqueda en cuadrícula para identificar el valor óptimo de CY que minimice alguna métrica de error del modelo (por ejemplo, MAPE, RMSE, etc.).

En el siguiente ejemplo, utilizaremos una búsqueda en cuadrícula para ajustar los parámetros del modelo oz y lo, que minimizan el MAPE del modelo para la serie de precios de Robusta. Por lo tanto, dividiremos el modelo en particiones de entrenamiento, prueba y validación, y para cada valor de a, aplicaremos los siguientes pasos:

1. Entrenar el modelo SES usando 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 prueba
3. Agregue las particiones de entrenamiento y 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 usar el conjunto de validación para validar el rendimiento del modelo. Antes de configurar la función, configuremos las particiones de entrenamiento, prueba y validación usando la función t s_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

Usaremos las variables trainl y testl para entrenar y probar la partición, 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 O . O 1 usando la función seq:

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

Dado que el valor de alfa debe ser mayor que cero, reemplazaremos 0 con un número pequeño que está bastante cerca de cero:

alpha[1] <- 0.001

Usaremos la función lapply para iterar en el modelo usando los diferentes valores de a:

library(dplyr)

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 puede ver en los siguientes resultados de prueba, mientras que a : 1 minimiza el MAPE en la partición de entrenamiento, oz : 0.03 minimiza la tasa de error en la partición de prueba. Los resultados de 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 = "purple", dash = "dot"), name = "Validation") %>%
  layout(title = "SES Model Grid Search Results",
         yaxis = list(title = "MAPE (%)"))

Figura 7. Como pudimos notar, esta grafica esta hecha con base a todo lo explicado anteriormente.

Metodo holt

El método de Holt, también conocido como modelo de suavizado exponencial doble, es una versión ampliada del modelo SES. Se basa en estimar el nivel y la tendencia más recientes con el uso de dos parámetros de suavizado, a y 5. Una vez que el modelo estima el nivel y la tendencia más recientes (LT y TT, respectivamente), los utiliza para construir el pronóstico de la serie.

Predicción con la función Holt

El paquete de pronóstico proporciona una implementación del modelo Holt con la función ho lt. Esta función inicializa automáticamente los valores de L1 y T1, e identifica los valores de a y 5 que minimizan el SSE del modelo en el conjunto de entrenamiento. En el siguiente ejemplo, recuperaremos los datos trimestrales del Producto Interno 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

## The . series is a ts object with 1 variable and 37 observations
## Frequency: 4
## Start time: 2010 1
## End time: 2019 1

ts_plot(gdp, 
        title = "Gross Domestic Product",
        Ytitle = "Billions of Dollars",
        Xtitle = "Source: U.S. Bureau of Economic Analysis / fred.stlouisfed.org")

Figura 8. Notará en la siguiente gráfica que la serie del PIB tiene una fuerte tendencia lineal y ningún componente estacional (ya que la serie está ajustada estacionalmente):

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 de precisión:

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 la tendencia de crecimiento y el pronóstico se pueden observar con la función test_forecast:

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

Figura 9. Si bien el modelo de Holt se diseñó para manejar series de tiempo con tendencia lineal, el argumento exponencial de la función de Holt brinda la opción de manejar series con tendencias exponenciales o decrecientes cuando se establece en TRUE. Para 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 estableceremos 5 en 0,75 (un enfoque más sólido para identificar el valor óptimo de 5 sería utilizar una búsqueda en 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

podemos graficar los valores ajustados y pronosticados contra los datos reales con la función test_forecastfuncfi0n$

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

Figura 10. Como puede ver, la tasa de error del segundo modelo de Holt es más equilibrada, donde la relación entre el error en el conjunto de prueba y entrenamiento es de 2,5 y 2,1 para las métricas RMSE y MAPE, respectivamente.

Modelo de Holt-Winters

Cerraremos este capítulo con el tercer y más avanzado modelo entre la familia de modelos de pronóstico de suavización exponencial: el modelo de Holt-Winters. El modelo Holt-Winters (HW) es una versión extendida del modelo Holt y puede manejar datos de series de tiempo con componentes de tendencia y estacionales. Pronosticar el componente estacional requería un tercer parámetro y ecuación más suaves, además de los de nivel y tendencia.

Tanto los componentes de tendencia como los estacionales podrían tener una estructura aditiva o de multiplicidad, lo que agrega cierta complejidad al modelo, ya que existen múltiples combinaciones posibles: -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.

Como hicimos anteriormente, crearemos particiones de entrenamiento y prueba usando los últimos 12 meses de la serie para evaluar el desempeño del modelo:

USgas_par <- ts_split(USgas, 12)

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

A continuación, usaremos el modelo de HoltWinters para pronosticar los últimos 12 meses de la serie (o el conjunto de prueba):

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 es pronosticar los próximos 1 2 meses (o los valores del conjunto de prueba) y evaluar el rendimiento del modelo con las funciones precision y test_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. Graficar los valores ajustados y pronosticados nos proporciona una mejor comprensión del rendimiento del modelo. Como podemos ver en el gráfico anterior, el modelo HW está haciendo un buen trabajo al capturar los patrones estacionales de ambas series. Por otro lado, al modelo le falta 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 en cuadrícula de manera similar a como lo hicimos con el modelo SES. En este caso, hay tres parámetros para optimizar: 01, fl y 7. El paquete TSstudio proporciona una función de búsqueda de cuadrícula personalizada basada en el enfoque de backtesting para entrenar una función HoltWinters.

Por razones de eficiencia, comenzaremos con una búsqueda superficial que incluya incrementos más grandes en los valores de los parámetros. Esto nos ayudará a reducir el área de búsqueda y luego aplicar una búsqueda más profunda en esas áreas. La búsqueda superficial incluirá pruebas retrospectivas durante 6 períodos diferentes usando una secuencia entre 0 y 1 con un incremento de O . 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'.

El resultado de la siguiente cuadrícula proporciona cualquier combinación de a, 13 y 7 con la tasa de error en cada período de prueba y su media general. La tabla ordena la media general del modelo de la mejor combinación a la peor y plot_grid proporciona una vista intuitiva del rango óptimo de valores de cada parámetro mediante el uso de un diagrama de paracords. De forma predeterminada, la función resalta los 10 modelos principales:

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

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

Notará en la salida del gráfico de búsqueda que el rango óptimo de a varía entre 0 . 1 y 0 5, :6 está entre 0 y O . l, y 7 está entre 0 . 2 y 0 . 4. Esto nos ayudará a establecer el rango de la búsqueda de hiperparámetros para la búsqueda de cuadrícula más profunda al reducir el rango de búsqueda pero usando 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'.

Como puede ver en la siguiente gráfica, el rango de error del 10% de los mejores modelos ha disminuido con respecto al de la búsqueda superficial:

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

Este gráfico se puede ver en una vista 3D (al realizar una búsqueda de tres parámetros): plot_grid(deep_grid, type = “3D”, top = 250) La salida se muestra en la siguiente captura de pantalla:

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

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)

El último paso de este proceso es extraer los valores de los parámetros de suavizado óptimos del modelo de cuadrícula en función de la búsqueda, volver a entrenar el modelo HW y utilizarlo para pronosticar los valores futuros de la serie:

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)

Usaremos las funciones de precisión y test_forecast para revisar el rendimiento del modelo HW que ha sido optimizado por una búsqueda de cuadrícula:

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

La gráfica de los valores ajustados y pronosticados frente a los valores reales se puede ver de la siguiente manera:

test_forecast(actual = USgas, 
              forecast.obj = fc_hw_grid, 
              test = test)

Figura 16. Como puede ver en el resultado de las funciones precision y test_forecast, la optimización del modelo con búsqueda en cuadrícula proporciona, en este caso, un aumento en el rendimiento del modelo al reducir la puntuación MAPE de 6,95 % a 5,92 %.

Resumen

En este capítulo, presentamos el uso de un promedio ponderado de observaciones pasadas para datos de series de tiempo de pronóstico. Comenzamos con un enfoque de pronóstico simplista e ingenuo con la función de promedio móvil. Aunque esta función se limita a pronósticos a corto plazo y solo puede manejar series de tiempo sin componentes estacionales ni de tendencia, proporciona contexto para funciones de suavizado exponencial. La familia de modelos de pronóstico de suavizado exponencial se basa en el uso de diferentes parámetros de suavizado, es decir, a, ti y W, para modelar los componentes principales de los datos de series temporales: nivel, tendencia y estacional, respectivamente. Las principales ventajas de las funciones de suavizado exponencial son su simplicidad, su bajo costo para la computación y su modularidad, lo que les permite manejar diferentes tipos de datos de series temporales, como tendencias lineales y exponenciales y componentes estacionales.