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

## 
## 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 ###################################
## # We noticed you have dplyr installed. 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 enter or source() into this session won't     #
## # work correctly.                                                             #
## #                                                                             #
## # All package code is unaffected because it is protected by the R namespace   #
## # mechanism.                                                                  #
## #                                                                             #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## # You can 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.                                #
## ################################### WARNING ###################################

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

Esta funcion es utilizada para, con algunos simples pasos, como un modelo predictivo. primero recapitulemos la estructura del SMA para series de tiempos suavizadas:

Grafico

para rolling windows de 2 lados y:

Grafico para rolling windows de cola izquierda

•Yt es la t observación de la serie Y con T observaciones, donde 1<t<T

•Yt’ representa el valor suavizado de Yt

•m representa la longitud de la ventana móvil, donde m≤T

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:

Grafico

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

Grafico

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

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

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, podemos generalizar la ecuación anterior para el pronóstico de la observación T+n de la serie (donde n ≥1):

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.

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

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, ademas, las ecuaciones previas pueden ser expandidas de las 2da predicciones y siguientes.

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

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

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

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

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

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

Grafico

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:

Grafico

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:

Entrenar el modelo SES usando la partición de entrenamiento y pronosticar las observaciones de la partición de prueba
Evaluar el rendimiento del modelo en la partición de entrenamiento y prueba
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
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, el cálculo de la más reciente el nivel y la tendencia de la serie se pueden obtener con las siguientes ecuaciones:

Grafico

y para tendencias multiplicativas:

Grafico

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:

Grafico

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:

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

•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. Las siguientes ecuaciones describen el modelo HW para una serie con componente estacional:

Grafico

Grafico

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.

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.00001,1,0.1),
                                            beta = seq(0.00001,1,0.1),
                                            gamma = seq(0.00001,1,0.1)), 
                        parallel = TRUE,
                        n.cores = 4)

## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.

## Warning: Detected creation of a 'multiprocess' future, which are defunct in
## future (>= 1.32.0) [2023-03-06]. Instead, specify either 'multisession'
## (recommended) or 'multicore'. If still used, 'multiprocess' becomes the same as
## 'sequential'.

## Warning: Strategy 'multiprocess' is deprecated in future (>= 1.32.0)
## [2023-03-06]. Instead, explicitly specify either 'multisession' (recommended)
## or 'multicore'.

La 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   gamma        1        2        3        4         5
## 1  0.70001 0.00001 0.90001 2.621390 2.184884 4.338277 4.401605  8.634916
## 2  0.70001 0.00001 0.80001 2.582018 2.352198 4.543953 4.710060  9.453419
## 3  0.70001 0.10001 0.90001 3.218790 2.809555 4.906356 4.150055  9.431933
## 4  0.60001 0.10001 0.60001 4.077523 1.935952 4.374445 4.114311  9.710598
## 5  0.60001 0.00001 0.60001 2.766514 2.636861 4.747206 4.420087  8.803496
## 6  0.60001 0.00001 0.70001 3.612816 2.990076 4.722938 4.084422  8.023298
## 7  0.60001 0.00001 0.50001 2.581400 2.653893 5.054241 4.765640  9.763542
## 8  0.60001 0.20001 0.60001 4.399353 3.012592 5.440290 3.986024  9.271970
## 9  0.70001 0.00001 0.70001 2.774289 2.872573 5.009110 4.983488 10.284533
## 10 0.50001 0.00001 0.40001 2.715483 2.976848 5.165077 4.557873  9.298461
##           6     mean
## 1  5.976851 4.692987
## 2  5.425901 4.844592
## 3  4.722678 4.873228
## 4  5.054252 4.877847
## 5  6.386860 4.960171
## 6  6.936269 5.061637
## 7  5.813575 5.105382
## 8  5.039624 5.191642
## 9  5.372686 5.216113
## 10 6.618207 5.221991

Con la función plot_grid nos proporciona una visión intuitiva de los rangos optimos de los valores de cada parametro. El modelo por default resalta el top 10% de modelos.

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

Gracias a estos graficos podemos saber que el rango optimo de alpha es de 0.5 a 0.7, para beta es de 0 a 0.2 y para gamma es entre 0.4 a 0.9. Estos valores nos ayudan a encontrar el rango de los hiperparametros de busqueda para la red profunda de busqueda estrechando el rango de busqueda, pero usando una busqueda mas puntual.

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

La grafica tambien la podemos ver de una forma 3D.

plot_grid(deep_grid, type = "3D", top = 250)

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

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

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

md_hw_grid <- HoltWinters(train,
                          alpha = deep_grid$alpha,
                          beta = deep_grid$beta,
                          gamma = deep_grid$gamma)
fc_hw_grid <- forecast(md_hw_grid, h = 12)
accuracy(fc_hw_grid, test)

##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 4.106435 117.0347 90.37810  0.2004600 4.387451 0.7849046 0.1526260
## Test set     3.124700 113.5063 97.20942 -0.2685092 3.899686 0.8442324 0.1506711
##              Theil's U
## Training set        NA
## Test set     0.3687748

Para finalizar, sacamos los valores optimos de los parametros suavizados de la red de datos de la busqueda profunda, mantenemos en el modelo HW y las usamos para predecir valores futuros de la serie.

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

grafica de los valores ajustados y predecidos vs los valores actuales

como se puede ver al optimiazar el modelo con redes de busqueda nos da una elevación de la presentación del modelo al reducir el puntaje de MAPE de un 6.95% a un 5.92%

Capitulo 10 Informe 3

Luz Arrieta

2023-04-17

Forecasting with exponentiial smoothing models

librerias

Weighted moving average

Forecasting with exponential smoothing

Simple exponential smoothing model

Forecasting the Ses function

Model optimization with grid set

Holt method

Forescasting with the Holt function

Holts-Winters Models