Rpubs3

Capítulo 8: Estrategias de Pronóstico

El Capítulo 8 se centra en “Estrategias de Pronóstico,” destacando los métodos y conceptos comunes para el pronóstico de datos de series temporales. Aquí tienes un resumen de los puntos clave:

Introducción al Pronóstico: Esta sección aborda varios métodos de pronóstico para series temporales, resaltando la importancia de entender patrones subyacentes como tendencias, estacionalidad y ciclos. Técnicas de Pronóstico: Se describen diferentes técnicas de pronóstico, haciendo hincapié en sus casos de uso apropiados: Pronóstico Ingenuo: Usar observaciones anteriores para predecir valores futuros. Es simple, pero limitado en escenarios con tendencias o estacionalidad. Promedio Móvil: Suavizar datos promediándolos durante un período definido, útil para eliminar ruido y resaltar tendencias. Suavizado Exponencial: Una técnica de suavizado más sofisticada, que permite mayor flexibilidad para manejar tendencias y estacionalidad. Modelos ARIMA: Estos integran componentes autorregresivos (AR), de diferenciación (I) y de promedio móvil (MA) para capturar patrones complejos. Evaluación y Comparación de Modelos de Pronóstico: Se discute la importancia de evaluar y comparar diferentes modelos para determinar su precisión y confiabilidad. Se utilizan métricas como el Error Porcentual Absoluto Medio (MAPE) para medir el rendimiento del modelo. Análisis de Residuales: El capítulo enfatiza la necesidad de verificar los residuales (la diferencia entre valores predichos y valores reales) para asegurar que los modelos no tengan patrones inexplicados. Esto garantiza que los pronósticos sean confiables e imparciales. Intervalos de Confianza: Se discuten como un método para estimar el rango dentro del cual se espera que caigan las observaciones futuras. Esto es importante para entender la incertidumbre en los pronósticos. En resumen, el Capítulo 8 proporciona una visión general completa de las diversas estrategias de pronóstico utilizadas en el análisis de series temporales, enfocándose en la importancia de seleccionar el método adecuado, evaluar el rendimiento del modelo y considerar la incertidumbre en los pronósticos. Si deseas más información detallada sobre algún tema o concepto específico, puedo proporcionarte más detalles.

# Chapter 8 Code

# -------- Code Chank 1 --------
library(TSstudio)

## Warning: package 'TSstudio' was built under R version 4.3.3

data(USgas)

ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

# Using the series time index to set the start and end point of each partiton
train <- window(USgas, 
                start = time(USgas)[1], 
                end = time(USgas)[length(USgas) - 12])

test <- window(USgas, 
               start = time(USgas)[length(USgas) - 12 + 1], 
               end = time(USgas)[length(USgas)])

ts_info(train)

##  The train series is a ts object with 1 variable and 226 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2018 10

ts_info(test)

##  The test series is a ts object with 1 variable and 12 observations
##  Frequency: 12 
##  Start time: 2018 11 
##  End time: 2019 10

# -------- Code Chank 2 --------
# The sample.out argument set the size of the testing partition
# (and therefore the training partition)
USgas_partitions <- ts_split(USgas, sample.out = 12)

train <- USgas_partitions$train
test <- USgas_partitions$test

ts_info(train)

##  The train series is a ts object with 1 variable and 226 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2018 10

ts_info(test)

##  The test series is a ts object with 1 variable and 12 observations
##  Frequency: 12 
##  Start time: 2018 11 
##  End time: 2019 10

# -------- Code Chank 3 --------
library(forecast)

## Warning: package 'forecast' was built under R version 4.3.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

md <- auto.arima(train)
# -------- Code Chank 4 --------
checkresiduals(md)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.95, df = 18, p-value = 0.1263
## 
## Model df: 6.   Total lags used: 24

# -------- Code Chank 5 --------
fc <- forecast(md, h = 12)
# -------- Code Chank 6 --------
accuracy(fc, test)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  5.844136  97.81626 73.42657 0.1170672 3.522348 0.6376860
## Test set     37.847885 103.22848 81.46603 1.3107987 3.261643 0.7075062
##                      ACF1 Theil's U
## Training set -0.004183172        NA
## Test set     -0.046708926 0.3404092

# -------- Code Chank 7 --------
test_forecast(actual = USgas,
              forecast.obj = fc,
              test = test) 

# -------- Code Chank 8 --------
library(forecast)

naive_model <- naive(train, h  = 12)
test_forecast(actual = USgas,
              forecast.obj = naive_model,
              test = test)

accuracy(naive_model, test)

##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set  -1.028444 285.6607 228.5084 -0.9218463 10.97123 1.984522
## Test set     301.891667 499.6914 379.1417  9.6798015 13.28187 3.292723
##                   ACF1 Theil's U
## Training set 0.3761105        NA
## Test set     0.7002486  1.499679

# -------- Code Chank 9 --------
snaive_model <- snaive(train, h = 12)
test_forecast(actual = USgas,
              forecast.obj = snaive_model,
              test = test)

accuracy(snaive_model, test)

##                    ME     RMSE      MAE      MPE     MAPE     MASE       ACF1
## Training set 33.99953 148.7049 115.1453 1.379869 5.494048 1.000000  0.4859501
## Test set     96.45000 164.6967 135.8833 3.612060 5.220458 1.180103 -0.2120929
##              Theil's U
## Training set        NA
## Test set     0.4289964

# -------- Code Chank 10 --------
md_final <- auto.arima(USgas)

fc_final <- forecast(md_final, h = 12)
# -------- Code Chank 11 --------
plot_forecast(fc_final,
              title = "The US Natural Gas Consumption Forecast",
              Xtitle = "Year",
              Ytitle = "Billion Cubic Feet")

# -------- Code Chank 12 --------
fc_final2 <- forecast(md_final, 
                      h = 60, 
                      level = c(80, 90))

plot_forecast(fc_final2,
              title = "The US Natural Gas Consumption Forecast",
              Xtitle = "Year",
              Ytitle = "Billion Cubic Feet")

# -------- Code Chank 13 --------
fc_final3 <- forecast_sim(model = md_final,
                          h = 60, 
                          n = 500) 
# -------- Code Chank 14 --------
library(plotly)

## Warning: package 'plotly' was built under R version 4.3.3

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

fc_final3$plot %>% 
  layout(title = "US Natural Gas Consumption - Forecasting Simulation",
         yaxis = list(title = "Billion Cubic Feet"),
         xaxis = list(title = "Year"))

# -------- Code Chank 15 --------
set.seed(1234)
# -------- Code Chank 16 --------
# Cargar los paquetes necesarios
library(TSstudio)  # Para manipulación de series temporales
library(forecast)  # Para modelos de series temporales

# Cargar los datos
data(USgas)

# Información de la serie temporal
ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

# Dividir los datos en conjuntos de entrenamiento y prueba
train <- window(USgas, 
                start = time(USgas)[1], 
                end = time(USgas)[length(USgas) - 12])

test <- window(USgas, 
               start = time(USgas)[length(USgas) - 12 + 1], 
               end = time(USgas)[length(USgas)])

# Entrenar el modelo con auto.arima
arima_model <- auto.arima(train)

# Verificar los residuos del modelo
checkresiduals(arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.95, df = 18, p-value = 0.1263
## 
## Model df: 6.   Total lags used: 24

# Hacer predicciones con el modelo ARIMA
forecast_12 <- forecast(arima_model, h = 12)

# Evaluar la precisión del modelo con el conjunto de prueba
accuracy(forecast_12, test)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  5.844136  97.81626 73.42657 0.1170672 3.522348 0.6376860
## Test set     37.847885 103.22848 81.46603 1.3107987 3.261643 0.7075062
##                      ACF1 Theil's U
## Training set -0.004183172        NA
## Test set     -0.046708926 0.3404092

# Visualizar la serie temporal con el pronóstico
plot_forecast(forecast_12,
              title = "Pronóstico del Consumo de Gas Natural en EE. UU.",
              Xtitle = "Año",
              Ytitle = "Miles de millones de pies cúbicos")

# Ajustar y predecir con el modelo completo
final_arima_model <- auto.arima(USgas)
final_forecast_60 <- forecast(final_arima_model, h = 60)

# Visualizar el pronóstico extendido
plot_forecast(final_forecast_60,
              title = "Pronóstico Extendido del Consumo de Gas Natural en EE. UU.",
              Xtitle = "Año",
              Ytitle = "Miles de millones de pies cúbicos")

# -------- Code Chank 17 --------
library(TSstudio)  # Para manipulación de series temporales
library(forecast)  # Para modelos de series temporales
library(ggplot2)  # Para visualizaciones con ggplot
library(plotly)  # Para gráficos interactivos

# Cargar los datos
data(USgas)

# Información de la serie temporal
ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

# Dividir los datos en conjuntos de entrenamiento y prueba
train <- window(USgas, 
                start = time(USgas)[1], 
                end = time(USgas)[length(USgas) - 12])

test <- window(USgas, 
               start = time(USgas)[length(USgas) - 12 + 1], 
               end = time(USgas)[length(USgas)])

# Entrenar el modelo ARIMA en el conjunto de entrenamiento
arima_model <- auto.arima(train)

# Comprobar residuos para evaluar el modelo
checkresiduals(arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.95, df = 18, p-value = 0.1263
## 
## Model df: 6.   Total lags used: 24

# Predecir usando el modelo ARIMA
forecast_12 <- forecast(arima_model, h = 12)

# Evaluar la precisión del modelo usando el conjunto de prueba
accuracy(forecast_12, test)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  5.844136  97.81626 73.42657 0.1170672 3.522348 0.6376860
## Test set     37.847885 103.22848 81.46603 1.3107987 3.261643 0.7075062
##                      ACF1 Theil's U
## Training set -0.004183172        NA
## Test set     -0.046708926 0.3404092

# Graficar el pronóstico usando ggplot2 y plotly para mayor flexibilidad
forecast_plot <- autoplot(forecast_12) +
  ggtitle("Pronóstico del Consumo de Gas Natural en EE. UU.") +
  xlab("Año") +
  ylab("Miles de millones de pies cúbicos")

# Hacer el gráfico interactivo con plotly
interactive_plot <- ggplotly(forecast_plot)

## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomForecast() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

# Mostrar el gráfico interactivo
print(interactive_plot)

# Si se requiere un pronóstico más largo
final_arima_model <- auto.arima(USgas)
final_forecast_60 <- forecast(final_arima_model, h = 60)

# Visualizar el pronóstico extendido con plotly
forecast_plot_60 <- autoplot(final_forecast_60) +
  ggtitle("Pronóstico Extendido del Consumo de Gas Natural en EE. UU.") +
  xlab("Año") +
  ylab("Miles de millones de pies cúbicos")

interactive_plot_60 <- ggplotly(forecast_plot_60)

## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomForecast() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues

# Mostrar el gráfico interactivo extendido
print(interactive_plot_60)

Capítulo 9: Pronóstico con Regresión Lineal

El Capítulo 9, titulado “Pronóstico con Regresión Lineal”, explora cómo utilizar modelos de regresión para analizar series temporales y hacer pronósticos. Aquí están los puntos clave que se abordan:

Conceptos Básicos de Regresión Lineal: Se explica el fundamento de la regresión lineal, describiendo cómo se pueden utilizar las líneas de tendencia para identificar y predecir comportamientos en datos temporales. El capítulo muestra cómo una línea de regresión puede ayudar a proyectar el futuro basándose en patrones pasados. Construcción de un Modelo de Regresión: El capítulo guía al lector a través del proceso de construcción de un modelo de regresión lineal. Esto incluye la selección de variables independientes y dependientes, y cómo se puede utilizar la regresión para predecir el valor futuro de una serie temporal. Verificación del Modelo: Se discuten técnicas para validar la precisión y la robustez del modelo de regresión. Esto incluye métodos para comprobar la distribución de los residuales y asegurar que el modelo no muestra tendencias inesperadas. Evaluación de la Eficacia del Modelo: El capítulo presenta métricas para medir el rendimiento del modelo, como el Coeficiente de Determinación (R²) y el Error Cuadrático Medio (MSE). Estos indicadores ayudan a determinar cuán bien se ajusta el modelo a los datos y qué tan preciso puede ser para futuros pronósticos. Extensión a Modelos Avanzados: Finalmente, el capítulo muestra cómo se puede combinar la regresión lineal con otros enfoques para mejorar la precisión del pronóstico. Esto incluye la incorporación de regresión polinómica y la regresión con errores de ARIMA para abordar datos más complejos. En resumen, el Capítulo 9 ofrece una guía completa para usar la regresión lineal en el contexto de series temporales, proporcionando herramientas para crear, validar y mejorar modelos de pronóstico.

# Chapter 9 Code

# -------- Code Chank 1 --------
library(TSstudio)

data("USgas")

ts_plot(USgas,
        title = "US Monthly Natural Gas consumption",
        Ytitle = "Billion Cubic Feet",
        Xtitle = "Year")

# -------- Code Chank 2 --------
ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

# -------- Code Chank 3 --------
ts_decompose(USgas)

# -------- Code Chank 4 --------
USgas_df <- ts_to_prophet(USgas)
# -------- Code Chank 5 --------
head(USgas_df)

##           ds      y
## 1 2000-01-01 2510.5
## 2 2000-02-01 2330.7
## 3 2000-03-01 2050.6
## 4 2000-04-01 1783.3
## 5 2000-05-01 1632.9
## 6 2000-06-01 1513.1

# -------- Code Chank 6 --------
USgas_df$trend <- 1:nrow(USgas_df)
# -------- Code Chank 7 --------
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

USgas_df$seasonal <- month(USgas_df$ds, label = T)
# -------- Code Chank 8 --------
head(USgas_df)

##           ds      y trend seasonal
## 1 2000-01-01 2510.5     1      Jan
## 2 2000-02-01 2330.7     2      Feb
## 3 2000-03-01 2050.6     3      Mar
## 4 2000-04-01 1783.3     4      Apr
## 5 2000-05-01 1632.9     5      May
## 6 2000-06-01 1513.1     6      Jun

# -------- Code Chank 9 --------
h <- 12 # setting a testing partition length

train <- USgas_df[1:(nrow(USgas_df) - h), ]

test <- USgas_df[(nrow(USgas_df) - h + 1):nrow(USgas_df), ]
# -------- Code Chank 10 --------
md_trend <- lm(y ~ trend, data = train)

summary(md_trend)

## 
## Call:
## lm(formula = y ~ trend, data = train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -547.2 -307.4 -153.2  333.1 1052.6 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1751.0074    52.6435   33.26  < 2e-16 ***
## trend          2.4489     0.4021    6.09 4.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 394.4 on 224 degrees of freedom
## Multiple R-squared:  0.1421, Adjusted R-squared:  0.1382 
## F-statistic: 37.09 on 1 and 224 DF,  p-value: 4.861e-09

# -------- Code Chank 11 --------
train$yhat <- predict(md_trend, newdata = train)

test$yhat <- predict(md_trend, newdata = test)
# -------- Code Chank 12 --------
library(plotly)
plot_lm <- function(data, train, test, title = NULL){
  p <- plot_ly(data = data, 
               x = ~ ds, 
               y = ~ y, 
               type = "scatter",
               mode = "line",
               name = "Actual") %>%
    add_lines(x =  ~ train$ds,
              y = ~ train$yhat,
              line = list(color = "red"),
              name = "Fitted") %>%
    add_lines(x =  ~ test$ds,
              y = ~ test$yhat,
              line = list(color = "green", dash = "dot", width = 3),
              name = "Forecasted") %>%
    layout(title = title,
           xaxis = list(title = ""),
           yaxis = list(title = "Billion Cubic Feet"),
           legend = list(x = 0.05, y = 0.95))
  return(p)
}
# -------- Code Chank 13 --------
plot_lm(data = USgas_df, 
        train = train, 
        test = test,
        title = "Predicting the Trend Component of the Series")

# -------- Code Chank 14 --------
mape_trend <- c(mean(abs(train$y - train$yhat) / train$y),
                mean(abs(test$y - test$yhat) / test$y))

mape_trend

## [1] 0.1644088 0.1299951

# -------- Code Chank 15 --------
md_seasonal <- lm(y ~ seasonal, data = train)

summary(md_seasonal)

## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -608.98 -162.34  -50.77  148.40  566.89 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2030.88      14.43 140.747  < 2e-16 ***
## seasonal.L   -480.00      50.24  -9.554  < 2e-16 ***
## seasonal.Q   1103.83      50.17  22.000  < 2e-16 ***
## seasonal.C     72.42      50.05   1.447 0.149389    
## seasonal^4    174.60      50.07   3.487 0.000592 ***
## seasonal^5    288.01      50.13   5.745 3.13e-08 ***
## seasonal^6    -44.67      50.09  -0.892 0.373460    
## seasonal^7   -187.91      49.96  -3.762 0.000218 ***
## seasonal^8     84.95      49.84   1.705 0.089706 .  
## seasonal^9     46.16      49.78   0.927 0.354828    
## seasonal^10    77.55      49.76   1.559 0.120587    
## seasonal^11   -11.09      49.75  -0.223 0.823856    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 216.9 on 214 degrees of freedom
## Multiple R-squared:  0.7521, Adjusted R-squared:  0.7394 
## F-statistic: 59.04 on 11 and 214 DF,  p-value: < 2.2e-16

# -------- Code Chank 16 --------
train$yhat <- predict(md_seasonal, newdata = train)
test$yhat <- predict(md_seasonal, newdata = test)

plot_lm(data = USgas_df, 
        train = train, 
        test = test,
        title = "Predicting the Seasonal Component of the Series")

# -------- Code Chank 17 --------
mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y),
                   mean(abs(test$y - test$yhat) / test$y))

mape_seasonal

## [1] 0.08628973 0.21502100

# -------- Code Chank 18 --------
md1 <- lm(y ~ seasonal + trend, data = train)

summary(md1)

## 
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -514.73  -77.17  -17.70   85.80  336.95 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1733.7153    17.0794 101.509  < 2e-16 ***
## seasonal.L  -498.1709    29.6354 -16.810  < 2e-16 ***
## seasonal.Q  1115.2951    29.5872  37.695  < 2e-16 ***
## seasonal.C    78.9835    29.5103   2.676  0.00802 ** 
## seasonal^4   175.6505    29.5196   5.950 1.09e-08 ***
## seasonal^5   285.0192    29.5568   9.643  < 2e-16 ***
## seasonal^6   -49.3611    29.5319  -1.671  0.09610 .  
## seasonal^7  -192.3050    29.4540  -6.529 4.77e-10 ***
## seasonal^8    81.8787    29.3835   2.787  0.00581 ** 
## seasonal^9    44.4849    29.3480   1.516  0.13106    
## seasonal^10   76.8636    29.3372   2.620  0.00943 ** 
## seasonal^11  -11.2755    29.3353  -0.384  0.70109    
## trend          2.6182     0.1305  20.065  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared:  0.9142, Adjusted R-squared:  0.9094 
## F-statistic: 189.2 on 12 and 213 DF,  p-value: < 2.2e-16

# -------- Code Chank 19 --------
train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)


plot_lm(data = USgas_df, 
        train = train, 
        test = test,
        title = "Predicting the Seasonal Component of the Series")

# -------- Code Chank 20 --------
mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
              mean(abs(test$y - test$yhat) / test$y))
mape_md1

## [1] 0.04774945 0.09143290

# -------- Code Chank 21 --------
md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)

summary(md2)

## 
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.882e+03  2.199e+01  85.568  < 2e-16 ***
## seasonal.L  -4.917e+02  2.530e+01 -19.438  < 2e-16 ***
## seasonal.Q   1.121e+03  2.525e+01  44.381  < 2e-16 ***
## seasonal.C   8.247e+01  2.518e+01   3.275  0.00123 ** 
## seasonal^4   1.763e+02  2.519e+01   6.999 3.35e-11 ***
## seasonal^5   2.835e+02  2.522e+01  11.243  < 2e-16 ***
## seasonal^6  -5.175e+01  2.520e+01  -2.054  0.04123 *  
## seasonal^7  -1.946e+02  2.513e+01  -7.741 3.97e-13 ***
## seasonal^8   8.030e+01  2.507e+01   3.203  0.00157 ** 
## seasonal^9   4.362e+01  2.504e+01   1.742  0.08293 .  
## seasonal^10  7.651e+01  2.503e+01   3.057  0.00253 ** 
## seasonal^11 -1.137e+01  2.503e+01  -0.454  0.65005    
## trend       -1.270e+00  4.472e-01  -2.840  0.00494 ** 
## I(trend^2)   1.713e-02  1.908e-03   8.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9341 
## F-statistic: 246.1 on 13 and 212 DF,  p-value: < 2.2e-16

# -------- Code Chank 22 --------
train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)


plot_lm(data = USgas_df, 
        train = train, 
        test = test,
        title = "Predicting the Seasonal Component of the Series")

mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
              mean(abs(test$y - test$yhat) / test$y))

mape_md2

## [1] 0.03688770 0.04212618

# -------- Code Chank 23 --------
USgas_split <- ts_split(USgas, sample.out = h)

train.ts <- USgas_split$train

test.ts <- USgas_split$test
# -------- Code Chank 24 --------
library(forecast)

md3 <- tslm(train.ts ~ season + trend + I(trend^2))

summary(md3)

## 
## Call:
## tslm(formula = train.ts ~ season + trend + I(trend^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.635e+03  3.224e+01  81.738  < 2e-16 ***
## season2     -3.004e+02  3.540e+01  -8.487 3.69e-15 ***
## season3     -4.841e+02  3.540e+01 -13.676  < 2e-16 ***
## season4     -9.128e+02  3.540e+01 -25.787  < 2e-16 ***
## season5     -1.099e+03  3.540e+01 -31.033  < 2e-16 ***
## season6     -1.118e+03  3.540e+01 -31.588  < 2e-16 ***
## season7     -9.547e+02  3.540e+01 -26.968  < 2e-16 ***
## season8     -9.322e+02  3.541e+01 -26.329  < 2e-16 ***
## season9     -1.136e+03  3.541e+01 -32.070  < 2e-16 ***
## season10    -1.046e+03  3.541e+01 -29.532  < 2e-16 ***
## season11    -8.001e+02  3.590e+01 -22.286  < 2e-16 ***
## season12    -2.618e+02  3.590e+01  -7.293 5.95e-12 ***
## trend       -1.270e+00  4.472e-01  -2.840  0.00494 ** 
## I(trend^2)   1.713e-02  1.908e-03   8.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9341 
## F-statistic: 246.1 on 13 and 212 DF,  p-value: < 2.2e-16

# -------- Code Chank 25 --------
r <- which(USgas_df$ds == as.Date("2014-01-01"))
USgas_df$s_break <- ifelse(year(USgas_df$ds) >= 2010, 1, 0)
USgas_df$s_break[r] <- 1
md3 <- tslm(USgas ~ season + trend + I(trend^2) + s_break, data = USgas_df)
summary(md3)

## 
## Call:
## tslm(formula = USgas ~ season + trend + I(trend^2) + s_break, 
##     data = USgas_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -469.25  -50.68   -2.66   63.63  275.89 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.661e+03  3.200e+01  83.164  < 2e-16 ***
## season2     -3.054e+02  3.448e+01  -8.858 2.61e-16 ***
## season3     -4.849e+02  3.448e+01 -14.062  < 2e-16 ***
## season4     -9.272e+02  3.449e+01 -26.885  < 2e-16 ***
## season5     -1.108e+03  3.449e+01 -32.114  < 2e-16 ***
## season6     -1.127e+03  3.450e+01 -32.660  < 2e-16 ***
## season7     -9.568e+02  3.450e+01 -27.730  < 2e-16 ***
## season8     -9.340e+02  3.451e+01 -27.061  < 2e-16 ***
## season9     -1.138e+03  3.452e+01 -32.972  < 2e-16 ***
## season10    -1.040e+03  3.453e+01 -30.122  < 2e-16 ***
## season11    -7.896e+02  3.497e+01 -22.577  < 2e-16 ***
## season12    -2.649e+02  3.498e+01  -7.571 9.72e-13 ***
## trend       -1.928e+00  4.479e-01  -4.304 2.51e-05 ***
## I(trend^2)   1.862e-02  1.676e-03  11.113  < 2e-16 ***
## s_break      6.060e+01  2.836e+01   2.137   0.0337 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109 on 223 degrees of freedom
## Multiple R-squared:  0.9423, Adjusted R-squared:  0.9387 
## F-statistic: 260.3 on 14 and 223 DF,  p-value: < 2.2e-16

library(UKgrid)

## Warning: package 'UKgrid' was built under R version 4.3.3

UKdaily <- extract_grid(type = "data.frame",
                        columns = "ND",
                        aggregate = "daily")

head(UKdaily)

##    TIMESTAMP      ND
## 1 2005-04-01 1920069
## 2 2005-04-02 1674699
## 3 2005-04-03 1631352
## 4 2005-04-04 1916693
## 5 2005-04-05 1952082
## 6 2005-04-06 1964584

ts_plot(UKdaily,
        title = "The UK National Demand for Electricity",
        Ytitle = "MW",
        Xtitle = "Year")

# -------- Code Chank 26 --------
ts_heatmap(UKdaily[which(year(UKdaily$TIMESTAMP) >= 2016),],
           title = "UK the Daily National Grid Demand Heatmap")

# -------- Code Chank 27 --------
library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

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

UKdaily <- UKdaily %>%
  mutate(wday = wday(TIMESTAMP, label = TRUE),
         month = month(TIMESTAMP, label = TRUE),
         lag365 = dplyr::lag(ND, 365)) %>%
  filter(!is.na(lag365)) %>%
  arrange(TIMESTAMP)
str(UKdaily)

## 'data.frame':    4939 obs. of  5 variables:
##  $ TIMESTAMP: Date, format: "2006-04-01" "2006-04-02" ...
##  $ ND       : int  1718405 1691341 1960325 2023886 2026204 2008422 1981175 1770440 1749715 2012865 ...
##  $ wday     : Ord.factor w/ 7 levels "Sun"<"Mon"<"Tue"<..: 7 1 2 3 4 5 6 7 1 2 ...
##  $ month    : Ord.factor w/ 12 levels "Jan"<"Feb"<"Mar"<..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ lag365   : int  1920069 1674699 1631352 1916693 1952082 1964584 1990895 2003982 1811436 1684720 ...

# -------- Code Chank 28 --------
start_date <- min(UKdaily$TIMESTAMP)


UK_ts <- ts(UKdaily$ND, 
            start = c(year(start_date), yday(start_date)),
            frequency = 365)

# -------- Code Chank 29 --------
h <-  365
UKpartitions <- ts_split(UK_ts, sample.out = h)
train_ts <- UKpartitions$train
test_ts <- UKpartitions$test

train_df <- UKdaily[1:(nrow(UKdaily) - h), ]
test_df <- UKdaily[(nrow(UKdaily) - h + 1):nrow(UKdaily), ]
# -------- Code Chank 30 --------
md_tslm1 <- tslm(train_ts ~ season + trend)
fc_tslm1 <- forecast(md_tslm1, h = h)
test_forecast(actual = UK_ts,
              forecast.obj = fc_tslm1,
              test = test_ts)

accuracy(fc_tslm1, test_ts)

##                         ME     RMSE       MAE        MPE     MAPE      MASE
## Training set -4.781286e-12 121551.4 100439.64 -0.5721337 6.294111 0.8418677
## Test set     -1.740215e+04 123156.6  96785.27 -1.8735336 7.160573 0.8112374
##                   ACF1 Theil's U
## Training set 0.5277884        NA
## Test set     0.5106681  1.071899

# -------- Code Chank 31 --------
md_tslm2 <- tslm(train_ts ~ season + trend + wday, data = train_df)
fc_tslm2 <- forecast(md_tslm2, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
              forecast.obj = fc_tslm2,
              test = test_ts)

accuracy(fc_tslm2, test_ts)

##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  8.172823e-12 70245.98 52146.79 -0.1738708 3.167605 0.4370853
## Test set     -1.764563e+04 80711.71 65373.21 -1.3715505 4.682071 0.5479470
##                   ACF1 Theil's U
## Training set 0.7513664        NA
## Test set     0.6075598   0.68445

# -------- Code Chank 32 --------
md_tslm3 <- tslm(train_ts ~ season + trend + wday + month + lag365, data = train_df)
fc_tslm3 <- forecast(md_tslm3, h = h, newdata = test_df)
test_forecast(actual = UK_ts,
              forecast.obj = fc_tslm3,
              test = test_ts)

accuracy(fc_tslm3, test_ts) 

##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -9.836939e-13 69904.92 52006.75 -0.1717563 3.163385 0.4359116
## Test set     -1.754061e+04 81783.55 65957.66 -1.3613252 4.722083 0.5528457
##                   ACF1 Theil's U
## Training set 0.7500146        NA
## Test set     0.6094666 0.6925767

# -------- Code Chank 33 --------
summary(md_tslm3)$coefficients %>% tail(1)

##           Estimate Std. Error   t value     Pr(>|t|)
## lag365 -0.06038328 0.01545495 -3.907051 9.490321e-05

# -------- Code Chank 34 --------
anova(md_tslm3)

## Analysis of Variance Table
## 
## Response: train_ts
##             Df     Sum Sq    Mean Sq    F value    Pr(>F)    
## season     364 1.4279e+14 3.9227e+11    73.5348 < 2.2e-16 ***
## trend        1 7.2634e+13 7.2634e+13 13615.7078 < 2.2e-16 ***
## wday         6 4.5009e+13 7.5016e+12  1406.2214 < 2.2e-16 ***
## month       11 1.3721e+11 1.2473e+10     2.3382  0.007266 ** 
## lag365       1 8.1432e+10 8.1432e+10    15.2650  9.49e-05 ***
## Residuals 4190 2.2352e+13 5.3345e+09                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# -------- Code Chank 35 --------
final_md <- tslm(UK_ts ~ season + trend + wday + month + lag365, 
                 data = UKdaily)
# -------- Code Chank 36 --------
checkresiduals(final_md)

## 
##  Breusch-Godfrey test for serial correlation of order up to 730
## 
## data:  Residuals from Linear regression model
## LM test = 3301.1, df = 730, p-value < 2.2e-16

# -------- Code Chank 37 --------
UK_fc_df <- data.frame(date = seq.Date(from = max(UKdaily$TIMESTAMP) + days(1), 
                                       by = "day", 
                                       length.out = h))
# -------- Code Chank 38 --------
UK_fc_df$wday <- factor(lubridate::wday(UK_fc_df$date, label = TRUE), ordered = FALSE)

UK_fc_df$month <- factor(month(UK_fc_df$date, label = TRUE), ordered = FALSE)

UK_fc_df$lag365 <- tail(UKdaily$ND, h)
# -------- Code Chank 39 --------
UKgrid_fc <- forecast(final_md, h = h, newdata = UK_fc_df)
# -------- Code Chank 40 --------
plot_forecast(UKgrid_fc,
              title = "The UK National Demand for Electricity Forecast",
              Ytitle = "MW",
              Xtitle = "Year")

Capítulo 10: Modelos de suavizado exponencial para el pronóstico de series temporales

En el capítulo 10 del libro, se abordan los modelos de suavizado exponencial para el pronóstico de series temporales, también conocidos como modelos de Holt-Winters. Estos modelos son útiles para datos con tendencias y patrones de estacionalidad.

El capítulo comienza explicando el principio básico del suavizado exponencial: asignar más peso a las observaciones recientes y menos a las pasadas, lo que permite una mayor sensibilidad a los cambios recientes en los datos.

Luego, el capítulo introduce los tres tipos principales de modelos de suavizado exponencial:

Suavizado Exponencial Simple: Adecuado para series sin tendencia ni estacionalidad. Este modelo utiliza un solo parámetro de suavizado para predecir la siguiente observación basada en las anteriores. Suavizado Exponencial Doble (Modelo de Holt): Ideal para series con tendencia. Este modelo añade un segundo parámetro para capturar el cambio en la tendencia a lo largo del tiempo. Suavizado Exponencial Triple (Modelo de Holt-Winters): Diseñado para series con tendencia y estacionalidad. Este modelo incluye un tercer parámetro para abordar los patrones estacionales, permitiendo una mejor adaptación a series con fluctuaciones estacionales regulares. A lo largo del capítulo, se muestra cómo ajustar estos modelos a diferentes tipos de datos, eligiendo parámetros óptimos y validando su precisión. Se discuten técnicas para seleccionar el modelo más adecuado según el tipo de serie temporal y cómo se pueden ajustar los parámetros para mejorar la precisión del pronóstico.

El capítulo concluye con ejemplos prácticos que demuestran cómo aplicar estos modelos en entornos reales y qué precauciones se deben tener para evitar el sobreajuste o la infraestimación de la variabilidad en los datos. Con esto, se proporciona una comprensión completa de cómo los modelos de suavizado exponencial pueden ser una herramienta poderosa para el pronóstico de series temporales.

# Chapter 10 Code

# -------- Code Chank 1 --------
library(TSstudio)

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

# -------- Code Chank 2 --------
robusta <- Coffee_Prices[,1]
# -------- Code Chank 3 --------
library(plotly)

ts_plot(robusta,
        title = "The Robusta Coffee Monthly Prices",
        Ytitle = "Price in USD",
        Xtitle = "Year")

# -------- Code Chank 4 --------
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 than the length of the series")}
  
  if(m > nrow(df)){
    stop("The length of the rolling window must be shorter than the length of the series")}
  if(!is.null(w)){
    if(length(w) != m){
      stop("The weight argument is not aligned with the length of the rolling window")
    } else if(sum(w) !=1){
      stop("The sum of the average weight is different than 1")
    }
  }
  
  # Setting the average weigths
  if(is.null(w)){
    w <- rep(1/m, m)
  }
  
  # Setting the data frame 
  #-----------------------
  # Changing the Date object column name
  names(df)[1] <- "date" 
  # Setting the training and testing partition 
  # according to the forecast horizon
  df$type <- c(rep("train", nrow(df) - h), 
               rep("test", h)) 
  
  # Spreading the table by the partition type
  df1 <- df %>% spread(key = type, value = y)
  
  # Create the target variable
  df1$yhat <- df1$train
  
  
  # Simple moving average function
  for(i in (nrow(df1) - h + 1):nrow(df1)){
    r <- (i-m):(i-1) 
    df1$yhat[i] <- sum(df1$yhat[r] * w) 
  } 
  
  # dropping from the yhat variable the actual values
  # that were used for the rolling window
  df1$yhat <- ifelse(is.na(df1$test), NA, df1$yhat)
  
  df1$y <- ifelse(is.na(df1$test), df1$train, df1$test)
  
  return(df1)
}
# -------- Code Chank 5 --------
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)
# -------- Code Chank 6 --------
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."))

# -------- Code Chank 7 --------
data(USgas)

USgas_df <- ts_to_prophet(USgas)
# -------- Code Chank 8 --------
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))
# -------- Code Chank 9 --------
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"))

# -------- Code Chank 10 --------
alpha_df <- data.frame(index = seq(from = 1, to = 15, by = 1),
                       power = seq(from = 14, to = 0, by = -1))

alpha_df$alpha_0.01 <- 0.01 * (1 - 0.01) ^ alpha_df$power
alpha_df$alpha_0.2 <- 0.2 * (1 - 0.2) ^ alpha_df$power
alpha_df$alpha_0.4 <- 0.4 * (1 - 0.4) ^ alpha_df$power
alpha_df$alpha_0.6 <- 0.6 * (1 - 0.6) ^ alpha_df$power
alpha_df$alpha_0.8 <- 0.8 * (1 - 0.8) ^ alpha_df$power
alpha_df$alpha_1 <- 1 * (1 - 1) ^ alpha_df$power

plot_ly(data = alpha_df) %>%
  add_lines(x = ~ index, y = ~ alpha_0.01, name = "alpha = 0.01") %>%
  add_lines(x = ~ index, y = ~ alpha_0.2, name = "alpha = 0.2") %>%
  add_lines(x = ~ index, y = ~ alpha_0.4, name = "alpha = 0.4") %>%
  add_lines(x = ~ index, y = ~ alpha_0.6, name = "alpha = 0.6") %>%
  add_lines(x = ~ index, y = ~ alpha_0.8, name = "alpha = 0.8") %>%
  add_lines(x = ~ index, y = ~ alpha_1, name = "alpha = 1") %>%
  layout(title = "Decay Rate of the SES Weights",
         xaxis = list(title = "Index"),
         yaxis = list(title = "Weight"))

# -------- Code Chank 11 --------
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

# -------- Code Chank 12 --------
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)))

# -------- Code Chank 13 --------
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)))

# -------- Code Chank 14 --------
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
# -------- Code Chank 15 --------
alpha <- seq(from = 0, to = 1, by = 0.01)
# -------- Code Chank 16 --------
alpha[1] <- 0.001
# -------- Code Chank 17 --------
library(dplyr)

ses_grid <- lapply(alpha, function(i){
  md1 <- md_accuracy1 <- md2 <- md_accuracy2 <- results <-  NULL
  md1 <- ses(train1, h = 12, alpha = i, initial = "simple")
  md_accuracy1 <- accuracy(md1, test1)
  
  md2 <- ses(train2, h = 12, alpha = i, initial = "simple")
  md_accuracy2 <- accuracy(md2, valid)
  
  resutls <- data.frame(alpha = i, 
                        train = md_accuracy1[9], 
                        test = md_accuracy1[10], 
                        valid = md_accuracy2[10])
  
}) %>% bind_rows()

# -------- Code Chank 18 --------
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 (%)"))

# -------- Code Chank 19 --------
library(Quandl)

## Warning: package 'Quandl' was built under R version 4.3.3

## Loading required package: xts

## Warning: package 'xts' was built under R version 4.3.3

## Loading required package: zoo

## Warning: package 'zoo' was built under R version 4.3.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Attaching package: 'xts'

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

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

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

# -------- Code Chank 21 --------
gdp_par <- ts_split(gdp, sample.out = 8)

train <- gdp_par$train
test <- gdp_par$test
# -------- Code Chank 22 --------
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

# -------- Code Chank 23 --------
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

# -------- Code Chank 24 --------
test_forecast(gdp, forecast.obj = fc_holt, test = test)

# -------- Code Chank 25 --------
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

# -------- Code Chank 26 --------
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

# -------- Code Chank 27 --------
test_forecast(gdp, forecast.obj = fc_holt, test = test)

# -------- Code Chank 28 --------
data(USgas)

decompose(USgas) %>% plot()

# -------- Code Chank 29 --------
USgas_par <- ts_split(USgas, 12)

train <- USgas_par$train
test <- USgas_par$test
# -------- Code Chank 30 --------
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

# -------- Code Chank 31 --------
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)

Capítulo 11: Modelos ARIMA

El capítulo 11 está dedicado a los modelos ARIMA, una de las herramientas más utilizadas para el análisis y pronóstico de series temporales. El modelo ARIMA combina componentes autorregresivos (AR), de promedio móvil (MA), y de integración (I), para modelar datos que pueden tener tendencias, patrones estacionales y otros comportamientos complejos.

Primero, el capítulo explica cómo identificar el orden del modelo ARIMA usando gráficos de función de autocorrelación (ACF) y función de autocorrelación parcial (PACF). Esto ayuda a determinar cuántos términos AR y MA son necesarios para capturar la estructura de la serie temporal.

Luego, se introducen las variantes del modelo ARIMA para datos con estacionalidad, conocidos como modelos SARIMA, que añaden componentes para capturar patrones estacionales. Estos modelos son especialmente útiles cuando se trabaja con datos que muestran comportamientos cíclicos a intervalos regulares, como ventas mensuales o producción industrial.

El capítulo incluye varios ejemplos prácticos, mostrando cómo aplicar modelos ARIMA a datos reales, y resalta la importancia de verificar los supuestos del modelo, como la normalidad de los errores y la independencia de los residuos. Se presenta el uso de herramientas como la prueba de Ljung-Box para confirmar que los residuos son ruido blanco, indicando que el modelo está capturando correctamente la estructura de la serie temporal.

También se aborda la importancia del diagnóstico del modelo y cómo ajustar los parámetros para mejorar el rendimiento. Por último, el capítulo destaca el uso de herramientas como la función auto.arima, que automatiza el proceso de búsqueda de los mejores parámetros para un modelo ARIMA.

En resumen, el capítulo 11 ofrece una guía completa para trabajar con modelos ARIMA, desde la identificación y construcción del modelo hasta su validación y ajuste, proporcionando a los lectores las habilidades necesarias para abordar una amplia gama de problemas de pronóstico de series temporales.

# Chapter 11 Code

# -------- Code Chank 1 --------
set.seed(12345)

stationary_ts <- arima.sim(model = list(order = c(1,0,0),
                                        ar = 0.5 ),
                           n = 500)
# -------- Code Chank 2 --------
library(TSstudio)

ts_plot(stationary_ts, 
        title = "Stationary Time Series",
        Ytitle = "Value",
        Xtitle = "Index")

# -------- Code Chank 3 --------
set.seed(12345)

non_stationary_ts <- arima.sim(model = list(order = c(1,1,0),ar = 0.3 ),n = 500)

ts_plot(non_stationary_ts, 
        title = "Non-Stationary Time Series",
        Ytitle = "Value",
        Xtitle = "Index")

# -------- Code Chank 4 --------
data(AirPassengers)

ts_plot(AirPassengers,
        title = "Monthly Airline Passenger Numbers 1949-1960",
        Ytitle = "Thousands of Passengers",
        Xtitle = "Year")

# -------- Code Chank 5 --------
ts_plot(diff(AirPassengers, lag = 1),
        title = "AirPassengers Series - First Differencing",
        Xtitle = "Year",
        Ytitle = "Differencing of Thousands of Passengers")

# -------- Code Chank 6 --------
ts_plot(diff(diff(AirPassengers, lag = 1), 12),
        title = "AirPassengers Series - First and Seasonal Differencing",
        Xtitle = "Year",
        Ytitle = "Differencing of Thousands of Passengers")

# -------- Code Chank 7 --------
ts_plot(diff(log(AirPassengers), lag = 1),
        title = "AirPassengers Series - First Differencing with Log Transformation",
        Xtitle = "Year",
        Ytitle = "Differencing/Log of Thousands of Passengers")

# -------- Code Chank 8 --------
set.seed(12345)

library(plotly)

p1 <- plot_ly()
p2 <- plot_ly()

for(i in 1:20){
  rm <- NULL
  rw <- arima.sim(model = list(order = c(0, 1, 0)), n = 500)
  p1 <- p1 %>% add_lines(x = time(rw), y = as.numeric(rw))
  p2 <- p2 %>% add_lines(x = time(diff(rw)), y = as.numeric(diff(rw)))
}
# -------- Code Chank 9 --------
p1 %>% layout(title = "Simulate Random Walk",
              yaxis = list(title = "Value"),
              xaxis = list(title = "Index")) %>%
  hide_legend()

# -------- Code Chank 10 --------
p2 %>% layout(title = "Simulate Random Walk with First-Order Differencing",
              yaxis = list(title = "Value"),
              xaxis = list(title = "Index")) %>%
  hide_legend()

# -------- Code Chank 11 --------
set.seed(12345)

ar2 <- arima.sim(model = list(order = c(2,0,0),
                              ar = c(0.9, -0.3)),
                 n = 500) 
# -------- Code Chank 12 --------
ts_plot(ar2, 
        title = "Simulate AR(2) Series",
        Ytitle = "Value",
        Xtitle = "Index")

# -------- Code Chank 13 --------
md_ar <- ar(ar2)
md_ar

## 
## Call:
## ar(x = ar2)
## 
## Coefficients:
##       1        2  
##  0.8851  -0.2900  
## 
## Order selected 2  sigma^2 estimated as  1.049

# -------- Code Chank 14 --------
par(mfrow=c(1,2))
acf(ar2)
pacf(ar2)

# -------- Code Chank 15 --------
set.seed(12345)

ma2 <- arima.sim(model = list(order = c(0, 0, 2),
                              ma = c(0.5, -0.3)),
                 n = 500) 
# -------- Code Chank 16 --------
ts_plot(ma2, 
        title = "Simulate MA(2) Series",
        Ytitle = "Value",
        Xtitle = "Index")

# -------- Code Chank 17 --------
md_ma <- arima(ma2, order = c(0,0,2), method = "ML")
md_ma

## 
## Call:
## arima(x = ma2, order = c(0, 0, 2), method = "ML")
## 
## Coefficients:
##         ma1      ma2  intercept
##       0.530  -0.3454     0.0875
## s.e.  0.041   0.0406     0.0525
## 
## sigma^2 estimated as 0.9829:  log likelihood = -705.81,  aic = 1419.62

# -------- Code Chank 18 --------
par(mfrow=c(1,2))
acf(ma2)
pacf(ma2)

# -------- Code Chank 19 --------
set.seed(12345)

arma <- arima.sim(model = list(order(1,0,2), ar = c(0.7), ma = c(0.5,-0.3)), n = 500)
# -------- Code Chank 20 --------
ts_plot(arma, 
        title = "Simulate ARMA(1,2) Series",
        Ytitle = "Value",
        Xtitle = "Index")

# -------- Code Chank 21 --------
arma_md <- arima(arma, order = c(1,0,2))
arma_md

## 
## Call:
## arima(x = arma, order = c(1, 0, 2))
## 
## Coefficients:
##          ar1     ma1      ma2  intercept
##       0.7439  0.4785  -0.3954     0.2853
## s.e.  0.0657  0.0878   0.0849     0.1891
## 
## sigma^2 estimated as 1.01:  log likelihood = -713.18,  aic = 1436.36

# -------- Code Chank 22 --------
par(mfrow=c(1,2))
acf(arma)
pacf(arma)

# -------- Code Chank 23 --------
arima(arma, order = c(1, 0, 1), include.mean = FALSE)

## 
## Call:
## arima(x = arma, order = c(1, 0, 1), include.mean = FALSE)
## 
## Coefficients:
##          ar1     ma1
##       0.4144  0.8864
## s.e.  0.0432  0.0248
## 
## sigma^2 estimated as 1.051:  log likelihood = -723,  aic = 1452

# -------- Code Chank 24 --------
arima(arma, order = c(2, 0, 1), include.mean = FALSE)

## 
## Call:
## arima(x = arma, order = c(2, 0, 1), include.mean = FALSE)
## 
## Coefficients:
##          ar1     ar2     ma1
##       0.3136  0.1918  0.9227
## s.e.  0.0486  0.0484  0.0183
## 
## sigma^2 estimated as 1.019:  log likelihood = -715.26,  aic = 1438.52

# -------- Code Chank 25 --------
arima(arma, order = c(1, 0, 2), include.mean = FALSE)

## 
## Call:
## arima(x = arma, order = c(1, 0, 2), include.mean = FALSE)
## 
## Coefficients:
##          ar1     ma1      ma2
##       0.7602  0.4654  -0.4079
## s.e.  0.0626  0.0858   0.0832
## 
## sigma^2 estimated as 1.014:  log likelihood = -714.27,  aic = 1436.54

# -------- Code Chank 26 --------
arima(arma, order = c(2, 0, 2), include.mean = FALSE)

## 
## Call:
## arima(x = arma, order = c(2, 0, 2), include.mean = FALSE)
## 
## Coefficients:
##          ar1     ar2     ma1      ma2
##       0.7239  0.0194  0.4997  -0.3783
## s.e.  0.2458  0.1257  0.2427   0.2134
## 
## sigma^2 estimated as 1.014:  log likelihood = -714.26,  aic = 1438.51

# -------- Code Chank 27 --------
library(forecast)
best_arma <- arima(arma, order = c(1, 0, 2), include.mean = FALSE)

checkresiduals(best_arma)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2) with zero mean
## Q* = 5.3129, df = 7, p-value = 0.6218
## 
## Model df: 3.   Total lags used: 10

# -------- Code Chank 28 --------
ar_fc <- forecast(md_ar, h = 100)
# -------- Code Chank 29 --------
plot_forecast(ar_fc, 
              title = "Forecast AR(2) Model",
              Ytitle = "Value",
              Xtitle = "Year")

# -------- Code Chank 30 --------
data(Coffee_Prices)

robusta_price <- window(Coffee_Prices[,1], start = c(2000, 1))
# -------- Code Chank 31 --------
ts_plot(robusta_price,
        title = "The Robusta Coffee Monthly Prices",
        Ytitle = "Price in USD",
        Xtitle = "Year") 

# -------- Code Chank 32 --------
acf(robusta_price)
# -------- Code Chank 33 --------
robusta_price_d1 <- diff(robusta_price)
# -------- Code Chank 34 --------
par(mfrow=c(1,2))

acf(robusta_price_d1)
pacf(robusta_price_d1)

# -------- Code Chank 35 --------
robusta_md <- arima(robusta_price, order = c(1, 1, 0))
# -------- Code Chank 36 --------
summary(robusta_md)

## 
## Call:
## arima(x = robusta_price, order = c(1, 1, 0))
## 
## Coefficients:
##          ar1
##       0.2780
## s.e.  0.0647
## 
## sigma^2 estimated as 0.007142:  log likelihood = 231.38,  aic = -458.76
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE     MAPE     MASE
## Training set 0.002595604 0.08432096 0.06494772 0.08104715 4.254984 1.001542
##                     ACF1
## Training set 0.001526295

# -------- Code Chank 37 --------
checkresiduals(robusta_md)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 26.896, df = 23, p-value = 0.2604
## 
## Model df: 1.   Total lags used: 24

# -------- Code Chank 38 --------
data(USgas)
# -------- Code Chank 39 --------
ts_plot(USgas,
        title = "US Monthly Natural Gas consumption", 
        Ytitle = "Billion Cubic Feet",
        Xtitle = "Year")

# -------- Code Chank 40 --------
USgas_split <- ts_split(USgas, sample.out = 12)

train <- USgas_split$train
test <- USgas_split$test
# -------- Code Chank 41 --------
par(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)

# -------- Code Chank 42 --------
USgas_d12 <- diff(train, 12)

ts_plot(USgas_d12,
        title = "US Monthly Natural Gas consumption - First Seasonal Difference", 
        Ytitle = "Billion Cubic Feet (First Difference)",
        Xtitle = "Year")

# -------- Code Chank 43 --------
USgas_d12_1 <- diff(diff(USgas_d12, 1))

ts_plot(USgas_d12_1,
        title = "US Monthly Natural Gas consumption - First Seasonal and Non-Seasonal Differencing", 
        Ytitle = "Billion Cubic Feet (Difference)",
        Xtitle = "Year")

# -------- Code Chank 44 --------
par(mfrow=c(1,2))
acf(USgas_d12_1, lag.max = 60)
pacf(USgas_d12_1, lag.max = 60)

# -------- Code Chank 45 --------
p <- q <- P <- Q <- 0:2
# -------- Code Chank 46 --------
arima_grid <- expand.grid(p,q,P,Q)
names(arima_grid) <- c("p", "q", "P", "Q")
arima_grid$d <- 1
arima_grid$D <- 1
# -------- Code Chank 47 --------
arima_grid$k <- rowSums(arima_grid)
# -------- Code Chank 48 --------
library(dplyr)

arima_grid <- arima_grid %>% filter(k <= 7)


# -------- Code Chank 51 --------
USgas_best_md <- arima(train, order = c(1,1,1), seasonal = list(order = c(2,1,1)))
USgas_best_md

## 
## Call:
## arima(x = train, order = c(1, 1, 1), seasonal = list(order = c(2, 1, 1)))
## 
## Coefficients:
##          ar1      ma1    sar1     sar2     sma1
##       0.4247  -0.9180  0.0132  -0.2639  -0.7449
## s.e.  0.0770   0.0376  0.0894   0.0834   0.0753
## 
## sigma^2 estimated as 10160:  log likelihood = -1292.96,  aic = 2597.91

# -------- Code Chank 52 --------
USgas_test_fc <- forecast(USgas_best_md, h = 12)
# -------- Code Chank 53 --------
accuracy(USgas_test_fc, test)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  6.081099  97.85701 73.36854 0.1298714 3.517097 0.6371821
## Test set     42.211253 104.79281 83.09943 1.4913412 3.314280 0.7216918
##                      ACF1 Theil's U
## Training set  0.004565602        NA
## Test set     -0.049999868 0.3469228

# -------- Code Chank 54 --------
test_forecast(USgas, 
              forecast.obj = USgas_test_fc,
              test = test)

# -------- Code Chank 55 --------
final_md <- arima(USgas, order = c(1,1,1), seasonal = list(order = c(2,1,1)))
# -------- Code Chank 56 --------
checkresiduals(final_md)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(2,1,1)[12]
## Q* = 30.173, df = 19, p-value = 0.04964
## 
## Model df: 5.   Total lags used: 24

# -------- Code Chank 57 --------
USgas_fc <- forecast(final_md, h = 12)
# -------- Code Chank 58 --------
plot_forecast(USgas_fc,
              title = "US Natural Gas Consumption - Forecast",
              Ytitle = "Billion Cubic Feet",
              Xtitle = "Year")

# -------- Code Chank 59 --------
USgas_auto_md1 <- auto.arima(train)
USgas_auto_md1

## Series: train 
## ARIMA(2,1,1)(2,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1    sar1     sar2     sma1
##       0.4301  -0.0372  -0.9098  0.0117  -0.2673  -0.7431
## s.e.  0.0794   0.0741   0.0452  0.0887   0.0830   0.0751
## 
## sigma^2 = 10446:  log likelihood = -1292.83
## AIC=2599.67   AICc=2600.22   BIC=2623.2

# -------- Code Chank 60 --------
USgas_auto_md2 <- auto.arima(train, 
                             max.order = 5,
                             D = 1,
                             d = 1,
                             ic = "aic",
                             stepwise = FALSE, 
                             approximation = FALSE)
USgas_auto_md2

## Series: train 
## ARIMA(1,1,1)(2,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1     sar2     sma1
##       0.4247  -0.9180  0.0132  -0.2639  -0.7449
## s.e.  0.0770   0.0376  0.0894   0.0834   0.0753
## 
## sigma^2 = 10405:  log likelihood = -1292.96
## AIC=2597.91   AICc=2598.32   BIC=2618.08

# -------- Code Chank 61 --------
df <- ts_to_prophet(AirPassengers) %>% setNames(c("date", "y"))

df$lag12 <- dplyr::lag(df$y, n = 12)

library(lubridate)

df$month <- factor(month(df$date, label = TRUE), ordered = FALSE)

df$trend <- 1:nrow(df)
# -------- Code Chank 62 --------
par <- ts_split(ts.obj = AirPassengers, sample.out = 12)

train <- par$train

test <- par$test
# -------- Code Chank 63 --------
train_df <- df[1:(nrow(df) - 12), ]
test_df <- df[(nrow(df) - 12 + 1):nrow(df), ]
# -------- Code Chank 64 --------
md1 <- tslm(train ~ season + trend + lag12, data = train_df)
# -------- Code Chank 65 --------
checkresiduals(md1)

## 
##  Breusch-Godfrey test for serial correlation of order up to 24
## 
## data:  Residuals from Linear regression model
## LM test = 83.147, df = 24, p-value = 1.903e-08

# -------- Code Chank 66 --------
md2 <- auto.arima(train, 
                  xreg = cbind(model.matrix(~ month,train_df)[,-1], 
                               train_df$trend, 
                               train_df$lag12), 
                  seasonal = TRUE, 
                  stepwise = FALSE, 
                  approximation = FALSE)

# -------- Code Chank 69 --------
fc1 <- forecast(md1, newdata = test_df)

fc2 <- forecast(md2, xreg = cbind(model.matrix(~ month,test_df)[,-1], 
                                  test_df$trend, 
                                  test_df$lag12))
# -------- Code Chank 70 --------
accuracy(fc1, test)

##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -1.894492e-15 13.99206 11.47136 -0.1200023 4.472154 0.3541758
## Test set      1.079611e+01 20.82782 18.57391  1.9530865 3.828115 0.5734654
##                   ACF1 Theil's U
## Training set 0.7502578        NA
## Test set     0.1653119 0.4052175

accuracy(fc2, test)

##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set   0.4117148  8.353629  6.397538  0.2571543 2.549682 0.2100998
## Test set     -11.2123706 17.928195 13.353910 -2.4898843 2.924174 0.4385521
##                      ACF1 Theil's U
## Training set  0.005966712        NA
## Test set     -0.325044929 0.3923795

Capítulo 12: Aprendizaje automático para el pronóstico de series temporales

¡Bienvenidos a la presentación sobre el capítulo 12! Este capítulo trata sobre el uso del aprendizaje automático para el pronóstico de series temporales. A través de la librería h2o, exploraremos cómo modelos como Random Forest, Gradient Boosting Machine (GBM) y AutoML pueden ayudarnos a predecir tendencias futuras con precisión.

Comenzamos analizando las ventajas y los desafíos del aprendizaje automático en comparación con enfoques tradicionales. El aprendizaje automático puede ofrecer alta precisión, pero también puede ser más costoso y complicado. Veremos cómo equilibrar estos aspectos.

Después, profundizamos en cómo entrenar modelos de aprendizaje automático para pronosticar series temporales. El capítulo destaca la importancia de crear buenas características (feature engineering) para maximizar el rendimiento de los modelos. Usamos un conjunto de datos sobre ventas de vehículos en Estados Unidos para mostrar cómo estos modelos pueden identificar patrones clave.

El capítulo también aborda la optimización de modelos. Aprendemos sobre la búsqueda de cuadrículas (grid search), validación cruzada y otras técnicas para encontrar los mejores hiperparámetros. Además, se muestra cómo automatizar el proceso de selección de modelos con la función h2o.automl.

Para finalizar, veremos resultados prácticos del uso de estos modelos y cómo se comparan entre sí en términos de precisión y eficiencia. El capítulo 12 es un recurso esencial para aquellos interesados en aplicar el aprendizaje automático al pronóstico de series temporales y aprender técnicas avanzadas para mejorar la exactitud de los modelos.

# Chapter 12 Code

# -------- Code Chank 1 --------
library(TSstudio)

data(USVSales)
# -------- Code Chank 2 --------
ts_info(USVSales)

##  The USVSales series is a ts object with 1 variable and 528 observations
##  Frequency: 12 
##  Start time: 1976 1 
##  End time: 2019 12

# -------- Code Chank 3 --------
ts_plot(USVSales,
        title = "US Total Monthly Vehicle Sales",
        Ytitle = "Thousands of Units",
        Xtitle = "Year")

# -------- Code Chank 4 --------
ts_decompose(USVSales)

# -------- Code Chank 5 --------
USVSales_detrend <- USVSales - decompose(USVSales)$trend

ts_seasonal(USVSales_detrend, type = "box")

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

## Warning: Ignoring 1 observations

# -------- Code Chank 6 --------

# -------- Code Chank 7 --------
ts_lags(USVSales, lags = c(12, 24, 36))

# -------- Code Chank 1 --------
# Cargar los paquetes necesarios
library(TSstudio)  # Para manipulación de series temporales
library(forecast)  # Para modelos de series temporales
library(dplyr)
library(lubridate)
library(plotly)
library(h2o)  # Para modelos con 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 http://docs.h2o.ai
## 
## ----------------------------------------------------------------------

## 
## Attaching package: 'h2o'

## The following objects are masked from 'package:lubridate':
## 
##     day, hour, month, week, year

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

# Inicializar h2o
h2o.init(max_mem_size = "16G")

## 
## H2O is not running yet, starting it now...
## 
## Note:  In case of errors look at the following log files:
##     C:\Users\Ale\AppData\Local\Temp\RtmpiCDbJP/h2o_Ale_started_from_r.out
##     C:\Users\Ale\AppData\Local\Temp\RtmpiCDbJP/h2o_Ale_started_from_r.err
## 
## 
## Starting H2O JVM and connecting:  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         5 seconds 826 milliseconds 
##     H2O cluster timezone:       America/Bogota 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.26.0.1 
##     H2O cluster version age:    4 years, 9 months and 11 days !!! 
##     H2O cluster name:           H2O_started_from_R_Ale_dfv347 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   14.21 GB 
##     H2O cluster total cores:    8 
##     H2O cluster allowed cores:  8 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         Amazon S3, Algos, AutoML, Core V3, Core V4 
##     R Version:                  R version 4.3.2 (2023-10-31 ucrt)

## Warning in h2o.clusterInfo(): 
## Your H2O cluster version is too old (4 years, 9 months and 11 days)!
## Please download and install the latest version from http://h2o.ai/download/

# -------- Code Chank 2 --------
# Datos y partición de conjuntos de entrenamiento y prueba
data(USVSales)

df <- ts_to_prophet(window(USVSales, start = c(2010, 1)))
names(df) <- c("date", "y")

# -------- Code Chank 3 --------
# Visualizar el subconjunto de datos
ts_plot(df,
        title = "US Total Monthly Vehicle Sales (Subset)",
        Ytitle = "Thousands of Units",
        Xtitle = "Year")

# -------- Code Chank 4 --------
# Preparación de datos con dplyr y lubridate
df <- df %>% 
  mutate(
    month = factor(lubridate::month(date, label = TRUE), ordered = FALSE),
    lag12 = lag(y, n = 12),
    trend = 1:nrow(df),
    trend_sqr = trend ^ 2
  ) %>%
  filter(!is.na(lag12))

# -------- Code Chank 5 --------
# Partición para entrenamiento y prueba
h <- 12  # Tamaño de la ventana de prueba
train_df <- df[1:(nrow(df) - h), ]
test_df <- df[(nrow(df) - h + 1):nrow(df), ]

# -------- Code Chank 6 --------
# Crear conjunto para pronóstico
forecast_df <- data.frame(
  date = seq.Date(from = max(df$date) + lubridate::month(1), 
                  length.out = h, by = "month"),
  trend = seq(from = max(df$trend) + 1, length.out = h, by = 1),
  trend_sqr = seq(from = max(df$trend_sqr) + 1, length.out = h, by = 1)
)

forecast_df$month <- factor(lubridate::month(forecast_df$date, label = TRUE), ordered = FALSE) 
forecast_df$lag12 <- tail(df$y, 12)

# -------- Code Chank 7 --------
# Entrenamiento de modelos de regresión lineal
lr <- lm(y ~ month + lag12 + trend + trend_sqr, data = train_df)

# -------- Code Chank 8 --------
# Evaluar el modelo y obtener el resumen
summary(lr)

## 
## Call:
## lm(formula = y ~ month + lag12 + trend + trend_sqr, data = train_df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -146.625  -38.997    0.111   39.196  112.577 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 442.05021   59.89832   7.380 1.23e-10 ***
## monthFeb    112.73160   34.16141   3.300 0.001439 ** 
## monthMar    299.20932   54.24042   5.516 4.03e-07 ***
## monthApr    182.52406   42.53129   4.292 4.88e-05 ***
## monthMay    268.75603   51.28464   5.240 1.24e-06 ***
## monthJun    224.66897   44.26374   5.076 2.41e-06 ***
## monthJul    177.88564   42.21898   4.213 6.49e-05 ***
## monthAug    241.63260   47.00693   5.140 1.86e-06 ***
## monthSep    152.99058   37.04199   4.130 8.76e-05 ***
## monthOct    125.16484   35.04896   3.571 0.000601 ***
## monthNov    127.97288   34.18772   3.743 0.000338 ***
## monthDec    278.67994   51.09552   5.454 5.21e-07 ***
## lag12         0.33906    0.10738   3.158 0.002236 ** 
## trend         9.07747    2.00504   4.527 2.03e-05 ***
## trend_sqr    -0.05587    0.01221  -4.576 1.69e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 59.6 on 81 degrees of freedom
## Multiple R-squared:  0.9198, Adjusted R-squared:  0.9059 
## F-statistic: 66.36 on 14 and 81 DF,  p-value: < 2.2e-16

# -------- Code Chank 9 --------
# Predicciones con el modelo lineal y cálculo de MAPE
test_df$yhat <- predict(lr, newdata = test_df)

mape_lr <- mean(abs(test_df$y - test_df$yhat) / test_df$y)  # MAPE del modelo de regresión lineal
mape_lr  # Mostrar MAPE

## [1] 0.03594578

# -------- Code Chank 10 --------
# Convertir los datos a formato h2o para usar modelos de aprendizaje automático
train_h <- as.h2o(train_df)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

test_h <- as.h2o(test_df)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

forecast_h <- as.h2o(forecast_df)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

# -------- Code Chank 11 --------
# Definir variables de entrada y salida para modelos de h2o
x <- c("month", "lag12", "trend", "trend_sqr")
y <- "y"

# -------- Code Chank 12 --------
# Entrenamiento de un modelo Random Forest
rf_md <- h2o.randomForest(
  training_frame = train_h,
  nfolds = 5,
  x = x,
  y = y,
  ntrees = 500,
  stopping_rounds = 10,
  stopping_metric = "RMSE",
  stopping_tolerance = 0.0001,
  seed = 1234
)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

# -------- Code Chank 13 --------
# Verificar importancia de variables
h2o.varimp_plot(rf_md)

# -------- Code Chank 15 --------
# Entrenamiento de modelos con búsqueda de hiperparámetros
hyper_params_rf <- list(
  mtries = c(2, 3, 4),
  sample_rate = c(0.632, 0.8, 0.95),
  col_sample_rate_per_tree = c(0.5, 0.9, 1),
  max_depth = seq(1, 30, 3),
  min_rows = c(1, 2, 5, 10)
)

search_criteria_rf <- list(
  strategy = "RandomDiscrete",
  stopping_metric = "RMSE",
  stopping_tolerance = 0.0001,
  stopping_rounds = 10,
  max_runtime_secs = 60 * 20
)

rf2 <- h2o.grid(
  algorithm = "randomForest",
  search_criteria = search_criteria_rf,
  hyper_params = hyper_params_rf,
  x = x,
  y = y,
  training_frame = train_h,
  ntrees = 5000,
  nfolds = 5,
  grid_id = "rf_grid",
  seed = 1234
)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%

# -------- Code Chank 16 --------
# Obtener el mejor modelo del grid
rf2_grid_search <- h2o.getGrid(
  grid_id = "rf_grid",
  sort_by = "rmse",
  decreasing = FALSE
)

rf_grid_model <- h2o.getModel(rf2_grid_search@model_ids[[1]])