10. Forecasting with exponential smoothing models

Forecasting with exponentiial smoothing models

En este capitulo se verá como se pueden hacer estimaciones y series de tiempo teniendo en cuenta funciones y librerias que reducen el tiempo y entregan la misma información.

libraries

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

Forecasting with moving average model

hemos visto como usar las funciones para recopilar y hacer modelos de series de tiempo, ahora se vera colo usar las mismas funciones para realizar pronosticos

The simple moving average

esta funcion usada para series de tiempo puede ser programada con ciertos pasos para la realizacion de pronosticos.

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

como vemos, esta grafica es la serie de tiempo del precio del cafe tipo robusta en estados unidos desde 1960 hasta el 2019, toda la informacion esta contenida en data(Coffee_Prices) estos datos se introducen en una función relacionada con el paquete TSstudio y luego se grafica con la función ts_plot.

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 that 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 weigth 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 weights
  if (is.null(w)) {
    w <- rep(1/m, m)
  }
  ### setting the data frame ###
  # Changing theDate 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)
  }

Como vemos, con la ayuda del paquete tidyr se hizo una función SMA (simple moving average) para su posterior uso en gráficas de pronostico y series de tiempo.

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

Ésta gráfica se hizo con el propósito de comprobar la función realizada con anterioridad, podemos ver cómo extrae datos de data(Coffee_prices) y los guarda para posteriormente realizar la gráfica con la función plot_ly.

Como podemos observar, no solo se trata de una serie de tiempo (Línea contínua azul), tambien realiza pronósticos en base al promedio del café robusta en los estados unidos (líneas discontínuas de varios colores)

Weighted moving average

La WMA es una versión extendida de la SMA y útil cuando los datos se correlacionan con los pronósticos.

data(USgas)
USgas_df <- ts_to_prophet(USgas)
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))

library(plotly)
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"))

Aquí también se usó la función sma_forecast pero debido a la relación estrecha de el consumo de gas natural con el pronóstico de uso es preferible usar WMA para mostrar esa correlación entre su “peso” (líneas discontinuas) y la serie de tiempo real (línea continua)

Forecasting with exponential smoothing

Las funciones de pronóstico usando la forma “exponential smoothing” son muy parecidas a las SMA vistas anteriormente, su principal diferencia es que esta forma de pronosticar series de tiempo se basa principalmente en recolectar información en cada uno de los puntos anteriores y en base a esa información aproximar como será el futuro, mientras que la SMA solo una m datos anteriores y no infinitos puntos como esta función, por esta razón es la más usada para pronosticar. En esta sección veremos tres métodos de pronosticar con esta función, el “simple exponential smoothing model”, el método “Holt” y por último el método “Holt-Winters”

Simple exponential smoothing model

El simple exponential smoothing model (SES) como su nombre dice, crea modelos pronósticos como se explicó con anterioridad pero en específico trabaja haciendo suposiciones, la principal que hace es que en todo momento la serie de tiempo es constante y puede hacer el pronóstico, por lo tanto, para series de tiempo que dependan de una época del año o que varíen con frecuencia no es tan útil.

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, x = ~ TRUE, y = ~ TRUE) %>%
  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"))

## Warning: Can't display both discrete & non-discrete data on same axis

## Warning: 'scatter' objects don't have these attributes: 'Y'
## Valid attributes include:
## 'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'fillpattern', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'

En esta gráfica podemos ver de lleno como funciona el modelo SES, pues al principio extrae los valores de alpha usados para modelar la serie de tiempo, lo hace mediante una operación matemática que se puede ver con claridad en el chunk anterior, posterior a esto, se usa la función plot_ly para graficar los datos obtenidos con anterioridad en la función alpha_df.

Forecasting with the ses function

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

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

En la anterior gráfica se demostró como usar la función ses para hacer modelos predictivos, lo principal es escoger el tiempo para hacer la estimación, en este caso se usó el precio mensual del café robusta de los anteriores 12 meses, luego se hizo la parte predictiva y se hizo una función para que guardara la información, posterior a esto y con ayuda de la función test_forecast se realizó la gráfica con los datos contenidos en la anterior función. Como podemos ver la gráfica nos muestra el precio real ( línea azul) con el pronóstico (línea verde) y como se relacionan una con otra.

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 esta grafica se puede observar como para modelos planos se puede usar otra funcion llamada plot_forecast que nos sirve para modelos predictivos con intervalos de confianza asociados a dicha prediccion.

Model optimization with grid search

En los casos que utilizemos la funcion ses se puede reemplazar por los “grid search” que cumple la misma funcion de ses y ademas ayuda a minimizar los errores y por ende disminuye las SSE (sum of square errors)

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

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

podemos ver como se usa las funciones train1 y test1 para probar la prediccion y train2 para guardar el modelo predictivo y guardar los datos validados por la funcion valid.

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)
  
  results <- data.frame(alpha = i,
                        train = md_accuracy1[9],
                        test = md_accuracy1[10],
                        valid = md_accuracy2[10])
}) %>% bind_rows()

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

Esta grafica es la condensación de todo el trabajo hecho por los dos anteriores chunks de codigo con todas las funciones hechas para probar la predicción en base a los datos utilizados y todo el almacenamiento de esta información en funciones para usarlas luego por plot_ly para generar una grafica que compara el entrenamiento de la predicción, el test de esta y su posterior validación en diferentes colores y lineas para su mejor entendimiento.

Holt method

Este metodo es una version extendida del modelo SES, se basa en usar los ultimos datos y estimarlos usando dos variables, alpha y beta. Cuando completa esta parte, usa la estimación para construir predicciones.

Forecastting with holt function

Para hcer modelos predictivos con el modelo Holt usando la libreria forecast y la funcion holt, la cual inicializa las estimaciones de los ultimos datos y encuentra los alpha y beta que minimizan el SSE para al final poder ser usado para graficar la prediccion.

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 = "US Gross Domestic Product",
        Ytitle = "Billions of Dollars",
        Xtitle= "Source: U.S. Bureau of Economic Analysis / fred.stlouisfed.org")

En esta grafica se usaron GDP (gross domestic product) contenidos en la federal reserve economic data (FRED) API. notemos como se hicieron todas las variables dichas con anticipacion y luego se uso la funcion ts_plot para graficarlos.

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

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

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

En esta grafica se puede observar como en los anteriores dos chunks de codigo se extrajeron las variables alpha y beta con el proposito de crear una prediccion y luego se comparo con la funcion accuracy, la diferencia se grafico en la anterior diagrama.

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

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

Para el caso de esta grafica se hizo algo parecido, se busco los valores de alpha y beta y posterior se uso la funcion accuracy para determinar la diferencia entre ellos y poder graficarlos junto a la actual prediccion. La unica diferencia que presenta es que para estos datos se usaron datos mas precisos para mejorar la prediccion.

Holt-Winters model

Esta modelo es mas preciso que el modelo Holt pues puede manejar series de tiempo de diferentes formas, por lo tanto, puede tener multiples variables y hacerlo de mejor forma.

data(USgas)

decompose(USgas) %>% plot()

Para esta grafica se usaron los datos de data(USgas) y la funcion decompose para diagnosticar la estructura los componentes de la serie de tiempo.

USgas_par <- ts_split(USgas, 12)

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

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

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)

En esta grafica se ve todo el potencial de el modelo holt-winters en la forma de obtener, medificar y resscribir los datos de entrada para poder obtener una prediccion en base a los datos de entrada similar al modelo holt

shallow_grid <- TSstudio::ts_grid(ts.obj = train, 
                        model = "HoltWinters",
                        optim = "MAPE",
                        periods = 6, window_length = NULL,
                        window_space = 6,
                        window_test = 12,
                        hyper_params = list(alpha = seq(0.0001,1,0.1),
                                            beta = seq(0.0001,1,0.1),
                                            gamma = seq(0.0001,1,0.1)), 
                        parallel = TRUE,
                        n.cores = 8)

## 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'.

shallow_grid$grid_df[1:10,]

##     alpha   beta  gamma        1        2        3        4         5        6
## 1  0.7001 0.0001 0.9001 2.621532 2.182118 4.339228 4.402514  8.629448 5.963847
## 2  0.7001 0.0001 0.8001 2.582051 2.355368 4.545089 4.711022  9.448610 5.418789
## 3  0.7001 0.1001 0.9001 3.217019 2.815053 4.903808 4.151560  9.438374 4.719588
## 4  0.6001 0.1001 0.6001 4.073554 1.934255 4.373182 4.114618  9.712710 5.051915
## 5  0.6001 0.0001 0.6001 2.769566 2.631529 4.747894 4.420550  8.796590 6.374720
## 6  0.6001 0.0001 0.7001 3.617912 2.984865 4.725269 4.084648  8.015492 6.925019
## 7  0.6001 0.0001 0.5001 2.581533 2.653442 5.053519 4.766254  9.757433 5.801817
## 8  0.6001 0.2001 0.6001 4.400772 3.021856 5.436146 3.983245  9.278232 5.041236
## 9  0.7001 0.0001 0.7001 2.773422 2.875788 5.009255 4.984433 10.280172 5.366689
## 10 0.5001 0.0001 0.4001 2.717662 2.971851 5.165682 4.558245  9.290997 6.606041
##        mean
## 1  4.689781
## 2  4.843488
## 3  4.874234
## 4  4.876706
## 5  4.956808
## 6  5.058867
## 7  5.102333
## 8  5.193581
## 9  5.214960
## 10 5.218413

plot_grid(shallow_grid)

## Warning in base::is.null(grid.obj$parameters$hyper_params[[i]]) ||
## grid.obj$parameters$hyper_params[[i]] == : 'length(x) = 10 > 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) = 10 > 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) = 10 > 1' in coercion to
## 'logical(1)'

En esta grafica vemos como el modelo funciona para varios usos luego de realizar el modelo predictivo, para este caso y el siguiente al usar la formula parallel = true la grafica presenta datos en paralelo que demuestran la prediccion hecha por el modelo y condensado en la funcion shallow_grid

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

Para esta grafica ocurre algo similar, luego de condensar la informacion en la funcion deep_grid se realiza la grafica en la funcion plot_grid en forma paralela. La unica diferencia son los datos de entrada para alpha, beta y gamma, esta es la razon del cambio de vista de la grafica.

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

Para esta grafica se usan los mismos parametros que para el anterior, la unica diferencia es la presentacion de la informacion, pues en la anterior se mostrabba en forma paralela y con lineas mientras que en esta grafica se usa un modelo 3D para mejorar la visualizacion.

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)

Para esta grafica se usaron parte de los dos modelos anteriores, pues se uso la funcion deep_grid y se condenso la informacion como la primera grafica.

Summary

fin de la presentación