R Markdown

En este capitulo se expandira el uso de smoothing functions y se dara introduccion a sus aplicaciones de predicción. Se empezará explicando sobre moving average models y modelos smoothing exponenciales simples; y luego se agregaran mas y mas capas al modelo, al tiempo que el modelo maneja series de tiempos mas complejas.

##Forecasting with exponentiial smoothing models

##librerias

library (forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library (h2o)
## 
## ----------------------------------------------------------------------
## 
## 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)
library (plotly)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.2
## 
## 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.3.2
## Loading required package: xts
## Loading required package: zoo
## 
## 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

A continuación veremos un ejemplo usando SMA usando la base de datos Coffee-Price.

Esta base de datos es un objeto mts y contiene los precios mensuales de los cafes Robusta y Arabica entre los años 1960 y 2018. A continuación extraeremos los precios del cafe Robusta y procederemos a graficarlo usando la libreria plotly.

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

Ahora, crearemos una funcion basica SMA usando funciones del paquete tidyr:

library(tidyr)
sma_forecast <- function (df, h, m, w = NULL) { 
  # Error handling
  if (h > nrow(df)) {
    stop("The length of the forecast horizon must be shorter 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)
  }

Los argumentos de las funciones son:

• df: serie de entradas en un formato de datos de 2 columnas, la primera es un objeto fecha y la segunda son los valores reales de la serie.

• h: Horizonte de las predicciones. Ejemplo: la funcion establece unas ultimas h observaciones como conjunto de prueba, esto nos ayuda a evaluar el rendimiento de nuestro modelo.

• m: longitud de la ventana rodante.

• w: Peso de los promedios, usando, por defecto, los promedios aritmeticos.

La funcion sma_forecast tiene los siguientes componentes:

• 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.

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

• 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 24 meses de la serie Robusta utilizando una ventana móvil de 1, 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)

Usando el paquete plotly graficaremos los diferentes moving average functions:

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

Observando lo arrojado en la grafica podemos decir que:

• Si la longitud de la ventana rodante es corta: 1) el rango que se pronostica es cercano a observaciones recientes. 2) la rapidez con la que las predicciones convergen en un valor constante.

• Si la longitud de la ventana es mayor: 1) cuanto mas va a tardar que los valores de prediciion converjan en un valor constante. 2) pueden controlar mejor los choques y valores atipicos.

• Un modelo de prediccion SMA con ventana movil de longitud de 1 es igual a un modelo de predicción 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 temporales, sin tendencia ni patrones estacionales. Esto afecta principalmente al promedio aritmético que suaviza el patrón estacional y se vuelve plano a largo plazo.

##Weighted moving average.

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.

Grafico:

Aquí, YT+n es el valor n pronosticado de la serie en el tiempo T+n, y w es un escalar de tamaño m que define el peso de cada observación en la ventana móvil. El uso del promedio ponderado brinda más flexibilidad ya que puede manejar series con un patrón estacional. En el siguiente ejemplo, usaremos la función sma_forecast que creamos anteriormente para pronosticar los últimos 24 meses del conjunto de datos USgas. En este caso, utilizaremos el argumento w para establecer el peso promedio y, por lo tanto, transformar la función de SMA a WMA. Como hicimos anteriormente, primero transformaremos la serie en formato data.frame con el función to_to_prophet:

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

Se usaron 2 formas:

• El WMA para aplicar todos los ponderados en los seasonal lag.

• El WMA para aplicarle una ponderación a el primer lag un 0.2 y al ultimo seasonal lag uno de 0.8

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

Ambos modelos capturaron las ocilaciones de la serie de tiempo hasta cierto punto.

mientras el WMA capta los componentes de temporada de la serie de tiempo, no puede capturar los trends devido al efecto del promedio. Este modelo comienza a perder eficiencia una vez que cruza el horizonte de la longitud de la frecuencia, en este caso seria mas de un año para una serie mensual.

##Forecasting with exponential smoothing.

Entre los modelos tradicionales de pronóstico de series de tiempo, las funciones de suavizado exponencial son uno de los enfoques de pronóstico más populares. Se basa en pronosticar los valores futuros de la serie promediando las observaciones pasadas de la serie. La principal distinción entre el suavizado exponencial y el enfoque de promedio móvil es que el primero promedia todas las observaciones de la serie, a diferencia del segundo, un subconjunto de m observaciones.

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

•Modelo de suavizado exponencial simple

• Modelo Holt

• Modelo Holt Winters

##Simple exponential smoothing model.

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 a lo largo del tiempo y, por lo tanto, este modelo es adecuado para series sin componentes de tendencia ni estacionales.

El principal atributo del modelo SES es la técnica del promedio ponderado, que se basa en el decaimiento exponencial de los pesos de observación según su distancia cronológica desde los primeros valores pronosticados. La tasa de caída de los pesos de observación se establece mediante alpha, el parámetro de suavizado del modelo.

Además, la función SES es una función escalonada, donde el valor n pronosticado del modelo se convierte en la entrada del próximo pronóstico de las próximas observaciones, n+1.

##Grafico

Donde:

YT +1 es el valor pronosticado de la observación T+1 para una serie con n observaciones (por ejemplo, T = 1,…,n)

YT es la T observación de la serie

a es el parámetro de suavizado del modelo, donde 0 < a <1

Yr es el valor pronosticado de la observación T en el paso T-1

Dado que el modelo asume que el nivel del modelo no cambia con el tiempo.

##Grafico

Como se sabe el objetivo es calcularlos niveles de la serie a base de las bases de datos, resultando en una predicción plana.

##Grafico

Esta formula la podemos seguir expandiendo agregandole los valores de predicciónes anteriores.Se expandera hasta que llega al primer valor de las predicciones (cronologicamente) y2:

ademas cronologicamente, la primera prediccion es Y2.

##Grafico

al no haber datos que soporten a y1, se inician los datos con el ^Y1.Esto es tipicamente hecho asignando el valor actual de la primera observación o por estimaciónes.

##Grafico

otra forma de interpretar las predicciones del modelo SES se observan manipulando la ecuacion de prediccion.

##Grafico

Yt-^Yt denota el término de error del valor ajustado de la observación t. Por lo tanto, el valor pronosticado de la observación T+1 por el modelo SES no es más que el valor pronosticado de la observación anterior, ^Yt, y la proporción del término de error de pronóstico en función de a.

##Grafico

eT se define como el error de pronóstico para T observación de la serie (por ejemplo, la última observación de una serie con T observaciones),

##Grafico

Una forma más común para esta ecuación es la forma de componentes, que en el caso del modelo SES es la estimación de YT+1 e Yr para denotar la estimación del modelo del nivel de la serie en el tiempo T y T-1, respectivamente.

##Grafico

El parámetro de suavizado del modelo, a, define la tasa de caída del peso del modelo. Dado que a está más cerca de 1 (n=1), los pesos de las observaciones más recientes son más altos. Por otro lado, la tasa de decaimiento, en este caso, es más rápida.

##Grafico

el siguiente ejemplo demuestra el decaimiento del promedio de las observaciones en las 15 ultimas observaciones para valores de alpha 0.01 a 1:

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'

Se transforman las ecuaciones anteriores en componente de que describe el modelo por sus componentes. En el caso del modelo SES, tenemos un solo componente: el nivel del modelo. Como mencionamos anteriormente, el supuesto principal del modelo es que el nivel de la serie se mantiene igual a lo largo del tiempo.

##Forecasting the Ses function.

inicial: define el método para inicializar el valor de 1, que se puede calcular usando las primeras observaciones de la serie estableciendo el argumento en simple, o estimándolo con el modelo ets 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 s 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

Como se puede ver del modelo, la función ses pone como valor inicial 0.69, que es bastante cercano a el valor de la primera observación(Y1=0,696), y pone el parametro alpha en 0.9999, lo que tecnicamente es a las predicciones Naive. Ahora probemos el modelo.

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

La función ses esta usando el set de entrenamiento para identificar el nivel de la serie estimando el parametro alpha y el nivel inicial del modelo(y1). El valor predecido se ve cerca del ultimo valorobservado desde la serie alpha, que en este caso es cercana a 1. Ya que el propocido del SES es predecir el nivel de la serie, en la grafica no se ven registradas las pequeñas ocilaciones.

En estimaciones planas, el intervalo de confianza del modelo es importante ya que el nivel de incertidumbre es alto. Tambien será util evaluar si los valores predecidos estan en el intervalo de confianza del modelo.

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

Tal como se ve en la grafica de preceding forecast, el set testeado esta en el rango de intervalo de confianza de 80%

##Model optimization with grid set.

La función SES optimiza los parámetros del modelo para minimizar la suma de los errores cuadráticos en el conjunto de entrenamiento. Una alternativa es usar una búsqueda en red para identificar los valores óptimos de los parámetros que minimizan el error del modelo. En este caso, se utilizará una búsqueda en red para encontrar el valor óptimo de O que minimice el MAPE del modelo para la serie de precios de Robusta. Se dividirá el modelo en particiones de entrenamiento, prueba y validación y se aplicarán estos pasos para cada valor de a. El modelo SES tiene un pronóstico plano, lo que crea una brecha entre los valores ajustados y el valor pronosticado:

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 en la partición de entrenamiento y prueba.

3.Agregar las particiones de entrenamiento y prueba (en orden cronológico) y volver a entrenar el modelo en la nueva partición 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 permiten 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 configurar la función, configuremos las particiones de entrenamiento, prueba y validación usando la función superior:

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

se usan el entrenamiento 1 y teste 1 para entrenar y testear y el train 2 lo usamos para la validación del modelo. los siguientes valores de alpha definiran el rango de busqueda. Como alpha es maypr que cero, le asignamos un valor pequeño cercano a este.

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

como se observa en el testeo con un alpha=1 se minimiza el MAPE en el training partition. mientras que alpha=0.03 se minimiza el error en el entrenamiento.

el output es visible en el siguiente grafico.

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

##Holt method.

El método de Holt 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 B. 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 usando la siguiente ecuación:

• Para una serie con tendencia aditiva, YT+h = LT + hTT

• Para una serie con tendencia multiplicativa, YT+1 = LT x hTr

Aquí, las variables de esas ecuaciones son las siguientes: YT+h es el valor pronosticado del h valor pronosticado de una serie con T observaciones.

• LT es la estimación del nivel de la serie en el tiempo T.

• TT es la estimación del impacto marginal de la tendencia de la serie en el momento T.

• h es el número de pasos de pronóstico desde el tiempo T.

El cálculo del nivel y la tendencia de la serie por el modelo de Holt se basa en un promedio ponderado con el uso de dos parámetros de suavizado, a y B. Para una serie con una tendencia aditiva.

##Forescasting with the Holt function.

El paquete de predicción 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 a y B que minimizan el SSE del modelo en el conjunto de entrenamiento. En el siguiente ejemplo, recuperaremos los datos trimestrales del (PIB) de EE. UU. de la API de Datos Económicos de la (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

Se nota en el siguente grafico que GDP tiene una linea trend bastante marcada y no tiene componentes temporales.

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

Como hicimos antes, dejamos una gran parte de los datos como testeo del modelo y entrenamos el modelo con lo que reste de observaciones de la serie usando la función **holt*. Tambien se revisan los parametros del modelo.

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

Se compara el los resultados del modelo en el entrenamiento con los del testeo haviendo uso de 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

Como puede ver en el resultado, la relación entre la tasa de error en el conjunto de prueba y entrenamiento es más de 15 veces para el RMSE y 4,5 para el MAPE. Esta gran proporción en las métricas de error se deriva principalmente de las siguientes dos razones:

•Los valores ajustados del modelo en el conjunto de entrenamiento no están limitados por una línea lineal (a diferencia de la salida del pronóstico).

•El crecimiento de la tendencia en los últimos trimestres pasa de una tasa de crecimiento lineal a una tasa exponencial.

Los cambios en la tendencia de crecimiento y el pronóstico se pueden observar con la función de pronóstico de prueba:

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

Mientra el modelo Holt es para manejar series de tiempo con tendencias centrales, el argumento exponencial de holt fue creado es usado para manejar series de tiempo con crecimientos exponenciales o decrecimientos cuando set es verdad.

para el siguiente caso usaremos el argumento exponencial para controlar el patron de crecimiento.Usando un ponderado mayor para este caso (beta= o.75).

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

se observa la exactitud del modelo:

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

Graficamos.

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

Algo a destacar del modelo es que no pude predecir cambios atipicos de el mundo, no se ajusta a lo que pasó en realidad debido a que el factor de que las ventas e ingresos bajaran ya no fue tanto por un factor de mercado si no algo expontaneo global como fue la pandemia.

##Holts-Winters Models.

El modelo HW 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 Los diferentes tipos de componentes 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

Por lo tanto, antes de construir un modelo HW, se necesita identificar la estructura de la tendencia y los componentes estacionales.

primero usemos una funcion para descomponer la estructura y temporadas de los componentes de la serie.

data(USgas)

decompose(USgas) %>% plot()

Los trend y componentes de temporada tienen una estructura aditiva. Ahora crearemos un modelo de etrenamiento para los ultimos 12 meses y asi evaluar el rendimiento.

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

Ahora usaremos el holt-winter model para predecir el comportamiento en los futuros 12 meses.

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)

Al 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.

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.

Para ser eficientes se empieza la busqueda shallow que incluyan lorgos incrementos en los valores de parametros.Esto se reaaliza en 6 periodos distintos en un rango de entre 0.00001 a 1 con incrementos de 0.1

como el codigo fue escrito antes de que se obsoletara la fucion multiprocess, adicionamos la libreria dragon para que nos de ayuda.

shallow_grid <- ts_grid (train, 
                         model = "HoltWinters",
                         periods = 6,
                         window_space = 6,
                         window_test = 12,
                         hyper_params = list(alpha = seq (0.00001,1,0.1),
                                             beta = seq (0.00001,1,0.1)),
                                             parallel = TRUE,
                                             n.cores = 1)

La siguiente tabla de resultado nos muestra todas las convinaciones de los parametros alpha, gamma y beta junto a su respectivo error de cada valor y un ponderado. La tabla esta organizada deforma que de primero encntremos el mejor promedio, mientras que abajo estará el pero promedio.

shallow_grid$grid_df[1:10,]
##      alpha    beta        1        2        3        4         5        6
## 1  0.60001 0.00001 2.572601 2.749792 4.858440 4.493930  9.535202 6.098639
## 2  0.70001 0.00001 3.434907 3.516453 5.403897 5.003323 10.866599 5.371996
## 3  0.10001 0.00001 2.918477 4.712904 6.337821 4.463508  6.692979 9.319122
## 4  0.20001 0.00001 2.986169 4.839133 6.479420 4.688269  8.088628 8.576504
## 5  0.10001 0.10001 3.400355 4.213983 6.051691 4.990198  7.513887 9.669555
## 6  0.30001 0.00001 4.082823 5.143506 6.547409 4.670105  8.824514 7.813415
## 7  0.50001 0.00001 6.227754 5.508103 5.388359 4.180177  8.789908 6.998465
## 8  0.80001 0.00001 4.315910 4.218212 5.885633 5.552441 12.126531 5.164581
## 9  0.40001 0.00001 5.286817 5.435407 6.359496 4.385449  8.822639 7.377918
## 10 0.10001 0.20001 3.210780 4.164858 6.139847 5.255890  9.060409 9.902695
##        mean
## 1  5.051434
## 2  5.599529
## 3  5.740802
## 4  5.943021
## 5  5.973278
## 6  6.180295
## 7  6.182128
## 8  6.210551
## 9  6.277954
## 10 6.289080