Forecasting with ExponentialSmoothing Models

Se prevé una ampliación en el empleo de funciones de suavizado y la introducción de sus aplicaciones en la predicción. Este proceso se iniciará con la explicación de modelos de promedio móvil y modelos de suavizado exponencial simple. A medida que se avance, se añadirán capas adicionales al modelo, lo que permitirá su manejo en el contexto de series de tiempo más complejas.

# Librerías usadas 

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

La función a contiunación es utilizada para, con algunos simples pasos, como un modelo predictivo. primero recapitulemos la estructura del SMA para series de tiempos suavizadas:

\[ Y^t_t= \frac{Y_{t-\frac{m-1}{2}}+...+Y_t+...+Y_{t-\frac{m-1}{2}} }{m} \] Por otro lado, para el windows de 2 lados:

\[ Y^t_t= \frac{Y_{t-m+1}+...+Y_{t-1}+Y_t}{m} \] Donde:

Generalmente con cualquier tipo de ventana móvil, incluye la observación t de la serie. La conversión de una función SMA en un modelo de pronóstico se basa en promediar las m observaciones consecutivas pasadas de la serie. Por ejemplo, Y+1|t=T, que es el valor de la primera observación de pronóstico de la serie (es decir, la observación T+1) es igual al promedio de las últimas m observaciones: Donde, los siguientes términos se utilizan en la ecuación anterior:

\[Y_{t+1|t=T}=\frac{Y_{t-m+1}+...+Y_{t-1}+Y_t}{m} \]

(Yt-m+1)+…+(Yt-1 )+Yt son las últimas m observaciones de la serie Y con T observaciones, dado t=T

Yt+1|t=T representa el primer valor pronosticado de una serie con T observaciones (es decir, dado que t es igual a la última observación de la serie T)

Al igual que la ecuación anterior, m representa la longitud de la ventana móvil

A partir de las fórmulas anteriores de la función de suavizado (Yt’) y las funciones de pronóstico (Yt+1|t=T) que el valor suavizado de la última observación de la serie (observación T) es el valor pronosticado de observal en T+1. Además, solo el primer valor pronosticado (observación T+1) se construye promediando solo los valores reales de la serie. A medida que cambiamos la ventana móvil para pronosticar las siguientes observaciones de la serie, los valores reales se reemplazan con los valores pronosticados previamente. Por ejemplo, el cálculo del segundo valor pronosticado en línea, T+2, se define mediante la siguiente expresión:

\[Y_{t+2|t=T}=\frac{Y_{t-m+2}+...+Y_{t-1}+Y_t+Y_{t+1}}{m} \] este modelo es recomendado cuando la serie no tiene una estructura de patrones como trends y seasonal. en este caso es mejor asumir que los valores de prediccion estan relativamente cercas a la ultima observacion de la serie.

A continuación se verá 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 se crea una función básica SMA por medio 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 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:

La funcion sma_forecast tiene los siguientes componentes:

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 bien la función SMA es bastante simple de usar y económica en potencia de cómputo, tiene algunas limitaciones:

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. La función WMA se puede formalizar mediante la siguiente expresión:

\[ Y_{T+n}=w_1Y_{T+n-m}+...+w_mY_{T+n-1} \]

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

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

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:

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. Podemos formalizar esta relación usando las siguientes ecuaciones:

\[Y_{T+1}=\alpha Y_T+(1-\alpha)Y_T \]

Donde:

Dado que el modelo asume que el nivel del modelo no cambia con el tiempo, podemos generalizar la ecuación anterior para el pronóstico de la observación T+n de la serie (donde n ≥1):

\[Y_{T+n}= Y_{T+1}= \alpha Y_T + (1- \alpha) Y_T \]

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

para identificar el nivel de el modelo se puede obtener de la proxima ecuación.

\[Y_{T+1}= \alpha Y_t+(1-\alpha)(\alpha Y_{T-1}+(1-\alpha)Y{T-1}=\alpha Y_T + \alpha(1-\alpha)Y_{T-1}+(1 - \alpha)^2 Y_{T-1})\] 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, ademas, las ecuaciones previas pueden ser expandidas de las 2da predicciones y siguientes.

\[Y_2=\alpha Y_1+(1-\alpha)Y_1\]

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. Se puede simplificar y comprimir mas, asignandole a la ecuación ^Y2 a la mano derecha de la ecuación.

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

\[Y_{T+1}=\alpha Y_T + (1- \alpha ) Y_T = \alpha Y_T + Y_T - \alpha Y_T \]

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, o como sigue:

\[Y_{T+1}=Y_T+ e\alpha_T \]

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), o como sigue:

\[e_T=Y_T- Y_T\]

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. La ecuación anterior se puede reescribir así:

\[Y_{T+1}=l_t=\alpha Y_T + (1-\alpha)l_{T-1} \]

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. Un caso especial es un modelo SES con a = 1, que es equivalente a un modelo de pronóstico ingenuo:

\[Y_{T+1|\alpha=1}=Y_T \]

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 y, por lo tanto, podemos reescribir la ecuación del modelo usando las siguientes notaciones:

\[Y_{T+1}=l_t \] Aquí, lT define el nivel del modelo en el momento T, que es un promedio ponderado de YT y el nivel del período anterior, lT-1:

\[Y_{T+1}=l_T=\alpha Y_T + (1- \alpha)l_{T-1} \]

Forecasting the Ses function

nicial: 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.

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, el cálculo de la más reciente el nivel y la tendencia de la serie se pueden obtener con las siguientes ecuaciones:

\[L_T=\alpha Y_T+\alpha(1-\alpha)(L_{$-1}+T_{T-1})\]

\[ T_T= \beta(L_T-L_{T-1})+(1- \beta)T_{T-1}\]

Y para las tendencias multiplicativas:

\[L_T=\alpha Y_T+\alpha(1-\alpha)(L_{T-1} \times T_{T-1} \]

\[T_T = \beta (\frac{L_T}{L_{T-1}})+(1-\beta)T_{T-1} \]

El modelo de Holt estima los valores de LT y TT, el nivel y la tendencia de la serie en el tiempo T, utilizando un promedio ponderado de todas las observaciones de la serie. Este proceso recursivo comienza con el pronóstico de la segunda observación del modelo:

\[L_2= \alpha Y_1 + \alpha(1-\alpha)(L_1+T_1)\]

\[T_2=\beta(L_2-L_1)+(1-\beta)T_1\]

L1 y L2 son el pronóstico del nivel de la primera y segunda observación de la serie, respectivamente.

T1 y T2 son el pronóstico de tendencia de la primera y la segunda observación de la serie, respectivamente.

Dado que no podemos pronosticar el valor de nivel y tendencia para la primera observación de la serie, tendremos que aproximar L1 y T2. El pronóstico de las próximas observaciones en línea se crea recursivamente utilizando la estimación de nivel y tendencia, así como los valores reales de las observaciones anteriores.

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:

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

Se grafica el resultado

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 lo 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:

Por lo tanto, antes de construir un modelo HW, se necesita identificar la estructura de la tendencia y los componentes estacionales. Las siguientes ecuaciones describen el modelo HW para una serie con componente estacional:

\[Y_{T+1}=L_T+hT_T+S_{T+h-m}\] \[L_T=\alpha(Y_T-S_{T-m})+(1-\alpha)(L_{T-1}+T_{T-1}) \] \[T_T=\beta(L_t-L_t)+(1-\beta)+(1-\beta)T_{T-1} \] \[S_T=\Upsilon(Y_T - L_T)+(1- \Upsilon)S_{T-m} \]

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 12meses.

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

El modelo esta aprendiendo principalmente de gamma y alpha (nivel y actualizaciones de temporada). no esta aprendiendo de beta. ahora taratermos de hacer que prediga los siguientes 12 meses.

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.