Hands_on_time_chapters11_12
Valentina Londoño
2023-10-25
Forecasting with ARIMA Models: Capítulo 11
La Media Movil Integrada Autoregresiva, o ARIMA es un modelo de predicción basado en la autoregresión y la media movil.
El proceso estacionario
Para hacer uso de ARIMA es necesario conocer que es el proceso estacionario. Este se basa en las siguientes condiciones:
- La media y la varianza de la serie no cambian en el tiempo
- La estructura de correlación de la serie, junto con sus “lags” permanece igual a través del tiempo
Para los siguientes ejemplos se utilizará la función arima.sim.
set.seed(12345)
stationary_ts <- arima.sim(model = list(order = c(1,0,0),
ar = 0.5 ),
n = 500)
library(TSstudio)
ts_plot(stationary_ts,
title = "Stationary Time Series",
Ytitle = "Value",
Xtitle = "Index")Figura 1: Serie de tiempo estacionaria
A partir de esta figura se puede observar que la media dela serie de tiempo no cambia demasiado, permanece cerca de la linea de 0. Esto representa una serie estacionaria.
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")Figura 2: Serie de tiempo no estacionaria
Esta figura es una representación de una serie de tiempo con una tendencia cambiante a través del tiempo.
Para que una serie de tiempo sea no estacionaria tiene que cumplir con:
*Tener una tendencia dominante: La media de la serie cambia a través del tiempo.
- Ser una serie con un componente estacionario multiplicativo: la varianza de la serie es función de una oscilación estacional a través del tiempo.
data(AirPassengers)
ts_plot(AirPassengers,
title = "Monthly Airline Passenger Numbers 1949-1960",
Ytitle = "Thousands of Passengers",
Xtitle = "Year")Figura 3: Numero mensual de pasajeros en una aerolinea entre 1949 y 1960
Esta figura es un claro ejemplo de la violación de los dos principios de la serie estacionaria. En ella se puede ver claramente un componente estacionario y un cambio en la tendencia.
Transformación de una serie no estacionaria en una serie estacionaria
Algunos de los métodos paraa transformar una serie no-estacionaria en una estacionaria son los siguientes:
Método de direfenciación de series de tiempo
Método de transformación logarítmica
Método del proceso de caminata aleatoria (Random Walk)
Metodo de difenciación de series de tiempo
La manera más común de transformar las series de tiempo no estacionarias en estacionarias es diferenciando sus lags. El propósito de esto es remover la tendencia para estabilizar la serie de tiempo.
Para esto se utilizará la función diff, la cual se encarga de diferencias el input de la serie con un lag específico
ts_plot(diff(AirPassengers, lag = 1),
title = "AirPassengers Series - First Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Thousands of Passengers")Figura 4: Serie de pasajeros aereos: primera diferenciación
En esta figura se puede evidenciar cómo la media de la serie es relativamente constante a través del tiempo. A pesar de esto, aún falta un criterio por cumplir: la eliminación de la variación de la serie de tiempo. Para ello se utiliza la diferencia estacional.
ts_plot(diff(diff(AirPassengers, lag = 1), 12),
title = "AirPassengers Series - First and Seasonal Differencing",
Xtitle = "Year",
Ytitle = "Differencing of Thousands of Passengers")Figura 5: Serie de pasajeros aereos: Primera diferenciación y diferenciación estacional
En esta figura se nota cómo la serie de tiempo fue estabilizada.
Método de transformación logarítmica
Este método busca utilizar logaritmos para transformar la serie de tiempo. Cabe resaltar que no es una sustitución para la diferenciación sino un complemento a esta.
ts_plot(diff(log(AirPassengers), lag = 1),
title = "AirPassengers Series - First Differencing with Log Transformation",
Xtitle = "Year",
Ytitle = "Differencing/Log of Thousands of Passengers")Figura 6: Serie de pasajeros aéreos: primera diferenciación y transformación logarítmica
En esta figura se nota cómo rápidamente la serie de tiempo se estabiliza en comparación al enfoque anterior.
Método del proceso de caminata aleatoria (Random Walk)
Este proceso describe la aleatoriedad de la serie de tiempo.
set.seed(12345)
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutp1 <- 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)))
}
p1 %>% layout(title = "Simulate Random Walk",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()Figura 7: Simulación de Random Walk
En este ejemplo se demuestra una simulación de 20 caminos diferentes de 500 pasos, todos empezando en otiempo 0 y punto 0.
p2 %>% layout(title = "Simulate Random Walk with First-Order Differencing",
yaxis = list(title = "Value"),
xaxis = list(title = "Index")) %>%
hide_legend()Figura 8: Simulación de Random Walk con diferenciación de primer orden
En esta figura se muestra el ejemplo anterior sin embargo se le añade la difernciación utilizada previamente. La serie de tiempo posee diferentes componentes, todos con medias y tendencias similares a través del tiempo.
Proceso AR
Este proceso define el valor actual de la serie, solo si esta es estacionaria.
En el siguiente ejemplo se utilizará la función arima.sim para simular el proceso de estructuración y posteriormente para amoldarse al modelo AR.
set.seed(12345)
ar2 <- arima.sim(model = list(order = c(2,0,0),
ar = c(0.9, -0.3)),
n = 500)
ts_plot(ar2,
title = "Simulate AR(2) Series",
Ytitle = "Value",
Xtitle = "Index")Figura 9: Simulación de la serie AR
La función ar permite adaptar la serie de tiempo para añadir un modelo AR y predecir valores futuros.
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.049Identificando el proceso AR y sus características
Cuando se trabaja en series de tiempo reales, es necesario identificar la estructura de la serie antes de modelarla. Para la identificación es necesario:
Categorizar el tipo de proceso: Si es AR, MA, ARMA
Identificar el orden del proceso: cuántos valores pasados utiliza para la predicción
par(mfrow = c(1, 2))
acf(ar2)
pacf(ar2)
La gráfica ACF determina cómo se relacionan los dato en determinado periodo de tiempo en relación a sus precedentes. Mientras que la PACF muestra la correlación entre dos observaciones. El propósito de esta gráfica es indicar cuáles son los valores más utiles para predecir valores futuros. La zona entre las lineas azules son los valores no significativos.
A partir de esta figura se puede establecer que la serie ar es un proceso de AR de segundo orden, o, AR (2)
Proceso de la media movil
En algunos casos el modelo no puede identificar todos los patrones por lol que alguna información queda como residual. El propósito del proceso de la media movil es capturar todos esos paterones residuales por medio del modelamiento de la relación entre el termino de error actual y los terminos ded error pasados.
En el proximo ejemplo se simulará la serie MA(2) con la función arima.sim
set.seed(12345)
ma2 <- arima.sim(model = list(order = c(0, 0, 2),
ma = c(0.5, -0.3)),
n = 500)
ts_plot(ma2,
title = "Simulate MA(2) Series",
Ytitle = "Value",
Xtitle = "Index")Figura 10: Simulación de la serie MA(2)
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.62Identificando el proceso MA y sus características
par(mfrow = c(1,2))
acf(ma2)
pacf(ma2)
A partir de esta gráfica es posible concluir que la serie ma2 pertenece a un proceso MA (2), o, MA de segundo orden.
Modelo ARMA
En algunos casos combinar los modelos AR y MA permite manejar datos de tiempo complejos.
Para esta estructura se utiliza la función arima.sim como será mostrado próximamente.
set.seed(12345)
arma <- arima.sim(model = list(order(1,0,2), ar = c(0.7), ma = c(0.5,-0.3)), n = 500)
ts_plot(arma,
title = "Simulate ARMA(1,2) Series",
Ytitle = "Value",
Xtitle = "Index")Figura 11: Simulación de la serie ARMA
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.36Identificando el proceso ARMA
par(mfrow=c(1,2))
acf(arma)
pacf(arma)
Prediciendo con los modelos AR, MA y ARMA
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoobest_arma <- arima(arma, order = c(1, 0, 2), include.mean = FALSE)
ar_fc <- forecast(md_ar, h = 100)
plot_forecast(ar_fc,
title = "Forecast AR(2) Model",
Ytitle = "Value",
Xtitle = "Year")Figura 13:Predicciones usando modelo AR(2)
Modelo ARIMA
Una de la slimitaciones de los modelos AR, MA y ARMA es qaue no pueden manejar series no estacionarias por lo que es necesario transformarlas. Para el modelo ARIMA se utiliza el proceso de Integrción para difernciar las series con sus lags.
Identificando el proceso ARIMA
Para esta sección se utilizará rl modelo arima, para ello es necesario:
*Identificar el grado de diferenciación que se requiere para transferir la serie a un estado estacionario
*Identidicar el proceso ARMA
Identificación del grado de diferenciación
data(Coffee_Prices)
robusta_price <- window(Coffee_Prices[,1], start = c(2000, 1))
ts_plot(robusta_price,
title = "The Robusta Coffee Monthly Prices",
Ytitle = "Price in USD",
Xtitle = "Year") Figura 14: Precio mensual del café Robusta desde el año 2000
Como se observa, los precios del café tienen una tendencia ascendente por lo que está en un estado estacionario.
Por medio de la función acf se intentará buscar la relación de la serie con sus lags.
acf(robusta_price)
Figura 15: Relación de la serie del precio de café Robusta con sus lags
Como se puede ver, la correlación de la serie con sus lags va decayendo a través del tiempo de forma linear.
robusta_price_d1 <- diff(robusta_price)
par(mfrow=c(1,2))
acf(robusta_price_d1)
pacf(robusta_price_d1)
Figura 16: Diferenciación de la serie de café Robusta
robusta_md <- arima(robusta_price, order = c(1, 1, 0))
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.001526295checkresiduals(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: 24Figura 17: Residuos obtenidos de ARIMA
A partir de esta figura se puede establecer que los residuos obtenidos de este modelo son insignificantes por lo qe es posible ignorarlos.
Modelo ARIMA Estacional (SARIMA) Por medio del siguiente se evidenciará lo aprendido en el capitulo 11del libro.
data(USgas)
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Figura 18: Consumo mensual de gas en Estados Unidos
Como se observa, esta serie de tiempo posee un fuerte patrón estacionalcon una tendencia creciente por lo que no es estacionaria por o que se utilizará el método SARIMA.
USgas_split <- ts_split(USgas, sample.out = 12)
train <- USgas_split$train
test <- USgas_split$test
par(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)
Figura 19: Relación entre los lags estacionales y las funciones ACF y PACF.
A partir de esto es posible determinar que la serie tiene una fuerte correlación con los lags estacionale sy no estacionales por lo que es necesario identificar si está en un estado estacionario o no.
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")Figura 20: Diferenciación estacional del consumo mensual de gas natural en Estados Unidos.
Esta figura se muestra inestable a pesar de haber eliminado la tendencia por lo que se intentará hacedr una primera diferenciación.
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")Figura 21: Primera diferenciación estacional y no eatcional del consumo mensual de gas natural en Estados Unidos.
En esta figura se muestra una estabilización de la serie alrededor del valor 0 por lo que puede empezar a considerarse una función estacionaria.
par(mfrow=c(1,2))
acf(USgas_d12_1, lag.max = 60)
pacf(USgas_d12_1, lag.max = 60)
Utilizando el siguiente código se debe de obtener un output del modelo SARIMA para establecer los diferentes parametrod en una tabla.
p <- q <- P <- Q <- 0:2
arima_grid <- expand.grid(p,q,P,Q)
names(arima_grid) <- c("p", "q", "P", "Q")
arima_grid$d <- 1
arima_grid$D <- 1
arima_grid$k <- rowSums(arima_grid)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionarima_grid <- arima_grid %>% filter(k <= 7)
arima_search <- lapply (1: nrow(arima_grid), function(i){
md <- NULL
md <- arima (train, order = c(arima_grid$p[i], 1, arima_grid$q [i]), seasonal = list(order = c(arima_grid$P[i],1, arima_grid$Q[i])), method = "CSS")
results <- data.frame(p= arima_grid$p[i], d = 1, q = arima_grid$q[i],
P = arima_grid$P[i], D = 1, Q = arima_grid$Q[i],
AIC = md$aic)
}) %>%
bind_rows() %>%
arrange(AIC)
head(arima_search)## p d q P D Q AIC
## 1 0 1 0 0 1 0 NA
## 2 1 1 0 0 1 0 NA
## 3 2 1 0 0 1 0 NA
## 4 0 1 1 0 1 0 NA
## 5 1 1 1 0 1 0 NA
## 6 2 1 1 0 1 0 NAEl código anterior debe de arrojar una tabla, sin embargo se produjo un error por lo que no fue posible ejecutarla.
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.91USgas_test_fc <- forecast(USgas_best_md, h = 12)
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.3469228test_forecast(USgas,
forecast.obj = USgas_test_fc,
test = test)Figura 22: Gas en Estados Unidos, valores eales vs. predicciones.
Como se puede ver, el modelo SARIMA captura las variaciones estacionales y de tendencia pero se le dificultan las variaciones estacionales.
A continuación se generará la predicción final con el modelo seleccionado.
final_md <- arima(USgas, order = c(1,1,1), seasonal = list(order = c(2,1,1)))
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: 24Figura 24: Residuos del modelo ARIMA
A partir de esta figura se puede determinar que los residups de esta serie son insignificantes y que la serie está normalmente distribuida por lo que la predicción final quedaría:
USgas_fc <- forecast(final_md, h = 12)
plot_forecast(USgas_fc,
title = "US Natural Gas Consumption - Forecast",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")Figura 25: Predicciones del consumo de gas mensual en Estados Unidos.
Regresión linear con errores ARIMA
df <- ts_to_prophet(AirPassengers)
names(df) <- c ("date", "y")
df$lag12 <- dplyr::lag(df$y, n = 12)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondf$month <- factor(month(df$date, label = TRUE), ordered = FALSE)
df$trend <- 1:nrow(df)
par <- ts_split(ts.obj = AirPassengers, sample.out = 12)
train <- par$train
test <- par$test
train_df <- df[1:(nrow(df) - 12), ]
test_df <- df[(nrow(df) - 12 + 1):nrow(df), ]
md1 <- tslm(train ~ season + trend + lag12 , data =train_df)
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-08Modelando los residuos con el modelo ARIMA
md2 <- auto.arima(train,
xreg = cbind(model.matrix(~ month, train_df) [, -1],
train_df$trend,
train_df$lag12),
seasonal = TRUE,
stepwise = FALSE,
approximation = FALSE)
summary(md2)## Series: train
## Regression with ARIMA(2,0,0)(2,0,0)[12] errors
##
## Coefficients:
## ar1 ar2 sar1 sar2 monthFeb monthMar monthApr monthMay
## 0.5849 0.3056 -0.4421 -0.2063 -2.7523 0.8231 0.2066 2.8279
## s.e. 0.0897 0.0903 0.1050 0.1060 1.8417 2.3248 2.4162 2.5560
## monthJun monthJul monthAug monthSep monthOct monthNov monthDec
## 6.6057 11.2337 12.1909 3.8269 0.6350 -2.2723 -0.9918
## s.e. 3.2916 4.2324 4.1198 2.9243 2.3405 2.4211 1.9172
##
## 0.2726 1.0244
## s.e. 0.1367 0.0426
##
## sigma^2 = 72.82: log likelihood = -426.93
## AIC=889.86 AICc=896.63 BIC=940.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4117148 8.353629 6.397538 0.2571543 2.549682 0.2100998
## ACF1
## Training set 0.005966712checkresiduals(md2)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(2,0,0)(2,0,0)[12] errors
## Q* = 25.816, df = 20, p-value = 0.172
##
## Model df: 4. Total lags used: 24Forecasting with Machine Learning Models: Capítulo 12
Estructura de la serie
library(TSstudio)
data(USVSales)
ts_plot(USVSales,
title = "US Total Montly Vehicle Sales",
Ytitle = "Thousands of Units",
Xtitle = "Year")Figura 26: Serie de tiempo de la ventatotal de vehículos en Estados Unidos entre 1976 y 2018
La siguiente figura es un diagrama de serie de tiempo representando las ventas totales de carros en Estados Unidos. En esta serie se aprecia un patrón repetido cada 10 años el cual muestra un decenso en las ventas a comienzo de la década para luego aumentar hasta la mitad de la década y posteriormente descender. Este diagrama muestra un indicador en la economía de EEUU a través del tiempo.
Componentes de la serie
Utilizando la función ts_decompose se pueden ver los componentes de la figura 1. Esta función descompone los factores que forman la gráfica total.
ts_decompose(USVSales)Figura 27: Descomposición de la serie de tiempo de la venta total de vehículos vendidos en Estados Unidos entre 1976 y 2018.
En esta figura se observan 4 diferentes series de tiempo, cada una de 4 factores diferentes. La primera es la serie de tiempo observada, es la secuencia de datos en determinado momento de tiempo.
Luego de esta está la tendencia de esta serie de tiempo la cual es el comportamiento que esta tiene a largo plazo; en ella se muestran los años 2000 y 2015 como la cúspide en la venta de vehículos y 1983 y 2010 como los años con mayor descenso).
En tercer lugar está la tendencia estacional que es el comportamiento de la serie de tiempo en cada año. Esta se muestra como un valor repetitivo cuyo punto más alto es 100 y el menor es -200.
De último está la aleatoriedad , esta es el cambio que se produce en los valores observados debido a fenómenos. La explicación de estos mismos es de compleja expplicación por lo que se toman como valores aleatorios.
Al adicionar estos 4 componentes se obtiene la figura 1
Análisis de temporadas
Para realizar un análisis de temporadas en la serie de tiempo, es necesario substraer de la serie, hacer una descomposición y utilizar la función ts_seasonal. Luego de esto será posible usar diagramas de caja y bigote para determinar las variaciones estacionales.
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 observationsFigura 28: Diagrama estacional de la venta de carros por medio de caja y bigote.
A partir de esta figura es posibble comparar las variaciones estacionales para cada mes entre 1976 y 2018. Aquí son utilizados diagramas de caja y bigote para cada mes.
Según esta figura es posible observar como durante todos los años el mes de enero presenta las ventas más bajas de todo el año, seguido por noviembre y luego febrero. Del mismo moso, la mayor cantidad de ventas son presentadas en los meses de marzo, mayo y junio para ver un descenso en los meses siguientes hasta llegar a diciembre.
De todos los meses el que presenta una mayor variación en la cantidad de vehículos vendidos alrededor de los años es el mes de diciembre mientras que el mes más “estable” es abril.
Análisis de correlación
ts_cor(USVSales)Figura 29: Correlación de la venta de vehículos en Estados Unidos con la primera estación
En esta figura se observa un análisis macro de la relación que
ts_lags(USVSales, lags = c(12,24,36))Figura 30: Correlación de la venta de vehículos en Estados Unidos con las primeras 3 estaciones
En esta figura se compara la relación que tiene la serie de tiempo con cada una de las primeras 3 estaciones.
Análisis exploratorio
A partir de la serie de tiempo es posible identificar que esta serie tiene una tendencia ciclica y se evidencia una estacionalidad según los meses. A pesar de esto, esta tendencia tiene una variabilidad cada 10 años por lo que es necesario usar la última década para pronosticar las variaciones futuras.
Para esto, se utiliza la función ts_to_prophet la cual se encarga de transformar una serie de tiempo en un data.frame
df <- ts_to_prophet(window(USVSales, start = c(2010,1)))
names(df) <- c("date", "y")
head(df)## date y
## 1 2010-01-01 712.469
## 2 2010-02-01 793.362
## 3 2010-03-01 1083.953
## 4 2010-04-01 997.334
## 5 2010-05-01 1117.570
## 6 2010-06-01 1000.455ts_plot(df,
title = "US Total Monthly Vehicle Sales (Subset)",
Ytitle = "Thousands of Units",
Xtitle = "Year")Figura 31: Ventas totales de vehículos en Estados Unidos desde 2010
En esta figura se toman los datos de los vehículos vendidos en Estados Unidos desde eol año 2010.
Se hace uso de esta gráfica para poder delimitar la información que será usada para generar predicciones de los años próximos.
Feature engineering
En esta sección se utilizará la serie de tiempo para generar características que permitan informar al modelo. Para ello serán usadas las tendencias de la serie, el componente estacional y la correlación.
Para esto se utilizarán los paquetes dplyr y lubridate.
library(dplyr)
library(lubridate)
df <- df %>% mutate(month = factor(month(date, label = TRUE), ordered = FALSE),
lag12 = lag(y, n = 12)) %>%
filter(!is.na(lag12))
df$trend <- 1:nrow(df)
df$trend_sqr <- df$trend ^ 2
str(df)## 'data.frame': 108 obs. of 6 variables:
## $ date : Date, format: "2011-01-01" "2011-02-01" ...
## $ y : num 836 1007 1277 1174 1081 ...
## $ month : Factor w/ 12 levels "Jan","Feb","Mar",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ lag12 : num 712 793 1084 997 1118 ...
## $ trend : int 1 2 3 4 5 6 7 8 9 10 ...
## $ trend_sqr: num 1 4 9 16 25 36 49 64 81 100 ...Con el código anteriormente expuesto es posible transformar la información original para permitir que nuestro modelo predictivo funcione de una manera ideal. De todas las variables que permiten que se genere la gráfica inicial, en esta sección se filtran todos los dato hasta quedar con 6 variables que permitirán la predicción. Estas son:
La fecha, el número de vehículos vendidos, el mes, la autocorrelación de los meses, la tendencia (variación de la serie de tiempo sin incluir las irregularidades) y la tendencia al cuadrado.
Evaluación del modelo
Para comparar diferentes modelos va a ser necesario comparar los inputs usados previamente. Para ello va a ser necesario probarlos.
Al tener un horizonte de pronóstico por 12 meses, el conjujnto de datos que será usado es el último año de ventas.
h <- 12
train_df <- df[1:(nrow(df) - h), ]
test_df <- df[(nrow(df) - h + 1) : nrow(df), ]Además de añadir los datos de los últimos 12 meses, es necesario añadir los inputs de la predicción en un data.frame nuevo.
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))
forecast_df$trend_sqr <- forecast_df$trend ^ 2
forecast_df$month <- factor(month(forecast_df$date, label = TRUE), ordered = FALSE)
forecast_df$lag12 <- tail(df$y, 12)Puntos de referencia del modelo
Para evaluar el rendimiento del modelo se debe comparar con un modelo base. Para este modelo en específico se entrenará una serie de tiempo con un modelo regresional para evaluar su rendimiento .
lr <- lm(y ~ month + lag12 + trend + trend_sqr, data = train_df)
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) 542.93505 72.59490 7.479 7.91e-11 ***
## 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 7.73667 1.72415 4.487 2.36e-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-16test_df$yhat <- predict(lr, newdata = test_df)
mape_lr <- mean(abs(test_df$y - test_df$yhat) / test_df$y)
mape_lr## [1] 0.03594578El valor obtenido con este código es el Error Porcentual Absoluto Medio o MAPE por sus siglas en inglés. Este valor es fundamental para medir la precisión del pronóstico.
Grupo h2o
Para esta sección se utilizará el paquete h2o. Este pauqete permite implementar algoritmos a escala, haciendo posible el machine learning.
library(h2o)##
## ----------------------------------------------------------------------
##
## Your next step is to start H2O:
## > h2o.init()
##
## For H2O package documentation, ask for help:
## > ??h2o
##
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit https://docs.h2o.ai
##
## ----------------------------------------------------------------------##
## Attaching package: 'h2o'## The following objects are masked from 'package: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, trunch2o.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\W11\AppData\Local\Temp\RtmpAD0JKa\file4a583945bd7/h2o_W11_started_from_r.out
## C:\Users\W11\AppData\Local\Temp\RtmpAD0JKa\file4a583cb7a23/h2o_W11_started_from_r.err
##
##
## Starting H2O JVM and connecting: .. Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 6 seconds 440 milliseconds
## H2O cluster timezone: America/Bogota
## H2O data parsing timezone: UTC
## H2O cluster version: 3.42.0.2
## H2O cluster version age: 3 months and 7 days
## H2O cluster name: H2O_started_from_R_W11_dfv347
## H2O cluster total nodes: 1
## H2O cluster total memory: 14.22 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
## R Version: R version 4.3.1 (2023-06-16 ucrt)train_h <- as.h2o(train_df)##
|
| | 0%
|
|======================================================================| 100%test_h <- as.h2o(test_df)##
|
| | 0%
|
|======================================================================| 100%forecast_h <- as.h2o(forecast_df)##
|
| | 0%
|
|======================================================================| 100%x <- c("month", "lag12", "trend", "trend_sqr")
y <- "y"h2o.init permite determinar el tamaño de memoria que se utilizará. Luego de esto se añade la información que será usada en el grupo.
Entrenar un modelo ML
El paquete h2o brinda herramientas para poder entrenar y probar los modelos.En esta sección se utilizará el modelo CV o de Validación Cruzada;este método está basado en:
- Dividir el set de entrenamiento
- Entrenar el modelo K veces. Se debe dejar un folder afuera luego de cada entrenamiento para poder probarlos
- Se alinea el modelo con los parámetros establecidos
- El modelo final se prueba
Prediciendo con el modelo Random Forest
Luego de preparar los datos y crear nuevas características, se puede construir el primer modelo predictor con el algoritmo Random Forest. Está basado en el componente Random, en el cual el input de cada modelo de arbol está basado en una muestra aleatoria; y el componente Forest que es la colección de modelos de arbol.
rf_md <- h2o.randomForest(training_frame = train_h,
nfolds = 5,
x = x,
y = y,
ntrees = 500,
stopping_rounds = 10,
stopping_metric = "RMSE",
score_each_iteration = TRUE,
stopping_tolerance = 0.0001,
seed = 1234)##
|
| | 0%
|
|=============== | 21%
|
|======================================================================| 100%h2o.varimp_plot(rf_md)
Figura 32: Importancia de variables: DRF
En este diagrama de barras se muestra la relevancia de cada variable en una escala de 0 a 1. De aquí es posible ver que la variable más importante para le model es lag12, seguido de trend_sqr, el mes y luego la tendencia con un valor de 0.5.
rf_md@model$model_summary## Model Summary:
## number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1 41 41 31590 8
## max_depth mean_depth min_leaves max_leaves mean_leaves
## 1 12 10.04878 45 66 56.70732En esta tabla se evidencia un resumen de la información del modelo. Se puede observar que se utilizaron 41 arboles de los 500 establecidos, todo esto como resultado de los parámetros usados.
library(plotly)
tree_score <- rf_md@model$scoring_history$training_rmse
plot_ly(x = seq_along(tree_score), y = tree_score,
type = "scatter", mode = "liine") %>%
layout(title = "The Trained Model Score History",
yaxis = list(title = "RMSE"),
xaxis = list(title = "NUm. of trees"))## Warning: Ignoring 1 observationsFigura 33: Historial de puntaje entrenado: Modelo Random Forest
test_h$pred_rf <- h2o.predict(rf_md, test_h)##
|
| | 0%
|
|======================================================================| 100%test_1 <- as.data.frame(test_h)
mape_rf <- mean(abs(test_1$y - test_1$pred_rf) / test_1$y)
mape_rf## [1] 0.04647164En esta secuencia se midió el rendimiento de los modelos. Se utilizó la función h2o.prediction para predecir los valores correspondientes de la serie. Del mismo modo, setranfirieron los datos a un data.frame para posteriormente calcular el MAPE.
El valor arrojado fue el puntaje de error que posee el modelo con su configuración preestablecida, siendo este un 4.6%.
Para reducir el porcentaje de error que presenta este modelo es necesario establecer perametros de búsqueda los cuales son indicados por:
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.0),
max_depth = c(seq(1, 30, 3)))Los parametros usados son:
- mtries: define las columnas para seleccionarlas aleatoriamente en cada nodulo del arbol
- sample_rate: establece el muestreo de filas para cada arbol
- max_depth: especifica la profudidad máxima del arbol
- min_rows: establece el número minimo de observaciones por cada hoja.
Al haber una mayor cantidad de parámetros,las predicciones serán mucho más precisas. Sin embargo, para mayor eficiencia al ser una serie de tiempo, lo ideal es hacer una búsqueda aleatoria.
El siguiente código busca añadir parámetros para que el modelo se vuelva más específico. A raíz de un error en la ejecución del código, no es posible ver la tabla.
search_criteria_rf <- list( strategy = "RandomDiscrete",
stopping_metric = "rmse", stopping_tolerance = 0.0001, stopping_rounds =
10, max_runtime_secs = 60 * 20 )rf2 <-h2o:: 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%
|
|====== | 9%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%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]])
test_h$rf_grid <- h2o.predict(rf_grid_model, test_h)##
|
| | 0%
|
|======================================================================| 100%mape_rf2 <- mean (abs(test_1$y -test_1$rf_grid) / test_1$y)
mape_rf2## [1] NaNEn adición a lo anteriormente realizado, también es ideal realizar el siguiente código para contribuir a la optimización del modelo.
plot_ly(data = test_1) %>%
add_lines(x = ~ date, y = ~y, name = "Actual") %>%
add_lines(x = ~ date, y = ~ yhat, name = "Linear Regression", line = list(dash = "dot")) %>%
add_lines(x = ~ date, y = ~ pred_rf, name = "Random Forest", line = list(dash = "dash")) %>%
add_lines(x = ~ date, y = ~ rf_grid_model , name = "Random Forest (grid)", line = list(dash = "dash")) %>%
layout(title = "Total Vehicle Sales - Actual vs. Prediction (Random Forest)",
yaxis = list(title = "Thousands of Units"),
xaxis = list(title = "Month"))## Warning in class(x) <- unique.default(c("AsIs", oldClass(x))): Setting class(x)
## to multiple strings ("AsIs", "H2ORegressionModel", ...); result will no longer
## be an S4 objectFigura 34: Venta total de vehículos: datos actuales vs. predicción usando Random Forest
Pronóstico con el modelo GBM
El algoritmo GBM es otro modelo basado en arboles. Utiliza boosting para entrenar diferentes subconjuntos de datos y repite el entrenamiento de los subconjuntos del modelo que hayan obtenido una alta tasa de error. Esto le permite al modelo aprender de errores pasados y mejorar la predicción.
En el siguiente ejemplo se demuestra el uso de la función h2o.gbm con el mismo conjunto de datos usados previamente
gbm_md <- h2o.gbm(
training_frame = train_h,
nfolds = 5,
x = x,
y = y,
max_depth = 20,
distribution = "gaussian",
ntrees = 500,
learn_rate = 0.1,
score_each_iteration = TRUE)##
|
| | 0%
|
|===== | 6%
|
|========= | 13%
|
|============= | 18%
|
|================= | 24%
|
|=================== | 28%
|
|====================== | 31%
|
|======================== | 34%
|
|========================= | 36%
|
|=========================== | 39%
|
|============================= | 41%
|
|=============================== | 44%
|
|================================ | 46%
|
|================================== | 49%
|
|==================================== | 52%
|
|====================================== | 54%
|
|======================================== | 57%
|
|=========================================== | 61%
|
|============================================= | 64%
|
|============================================== | 66%
|
|================================================ | 68%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|====================================================== | 77%
|
|======================================================= | 79%
|
|======================================================== | 81%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 85%
|
|============================================================= | 87%
|
|============================================================== | 89%
|
|=============================================================== | 91%
|
|================================================================= | 93%
|
|=================================================================== | 95%
|
|==================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100% h2o.varimp_plot(gbm_md)
Figura 35: Importancia de variables de 0 a 1 según modelo GBM
En esta figura se puede ver cómo la variable de mayor importancia es lag12 con un valor de 1. A pesar de que en ambos modelos la variable más importante es esta, la diferencia entre ambos es que en gbm la fundamental es lag12, la importancia de las demás variables es muy mínima; mientras que en Random Forest el resto de variables poseían una relevancia media.
test_h$pred_gbm <- h2o.predict (gbm_md, test_h)##
|
| | 0%
|
|======================================================================| 100%test_1 <- as.data.frame(test_h)
mape_gbm <- mean(abs(test_1$y - test_1$pred_gbm) /test_1$y)
mape_gbm## [1] 0.03857898El siguiente código es para ver la comparación entre los valores actuales vs. la predicción con GBM en las ventas totales de vehículos.
plot_ly(data = test_1) %>%
add_lines(x = ~ date, y = ~y, name = "Actual") %>%
add_lines(x = date , y = ~yhat, name = "Linear Regression", line = list(dash = "dot")) %>%
add_lines(x = ~date, y = ~ pred_gbm, name = "Graident Boosting Machine", line = list(dash = "dash")) %>%
layout(title = "Total Vehicle Sales - Actual vs. Prediction (Gradient Boosting Machine)",
yaxis = list(title = "Thousan of Units"),
xaxis = list (title ="Month"))Figura x: Venta total de vehículos: actual vs predicción
##Prediciendo con el modelo AutoML
autoML1 <- h2o.automl(training_frame = train_h,
x = x,
y = y,
nfolds = 5,
max_runtime_secs = 60*20,
seed = 1234)##
|
| | 0%
|
|== | 3%
## 21:19:32.557: AutoML: XGBoost is not available; skipping it.
## 21:19:36.207: _min_rows param, The dataset size is too small to split for min_rows=100.0: must have at least 200.0 (weighted) rows, but have only 96.0.
|
|====== | 8%
|
|======= | 10%
|
|=========== | 16%
|
|============= | 18%
|
|================ | 22%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%test_h$pred_autoML <- h2o.predict(autoML1@leader, test_h)##
|
| | 0%
|
|======================================================================| 100%test_1 <- as.data.frame(test_h)
mape_autoML <- mean(abs(test_1$y - test_1$pred_autoML) / test_1$y)
mape_autoML## [1] 0.03860951plot_ly(data = test_1) %>%
add_lines(x = ~ date, y = ~y, name = "Actual") %>%
add_lines(x = ~ date, y = ~ yhat, name = "Linear Regression", line = list(dash = "dot")) %>%
add_lines(x = ~ date, y = ~ pred_autoML, name = "autoML", line = list(dash = "dash")) %>%
layout(title = "Total Vehicle Sales - Actual vs. Prediction (Auto ML Model)",
yaxis = list(title = "Thousands of Units"),
xaxis = list(title = "Month"))EL MODELO FINAL
Del mismo modo, este código no pudo ser ejecutado pero este propone una predicción final desdel el año 2019 hasta el 2020.
forecast_h$pred_gbm <- h2o.predict(gbm_md, forecast_h)##
|
| | 0%
|
|======================================================================| 100%forecast_h$pred_rf <- h2o.predict(rf_grid_model, forecast_h)##
|
| | 0%
|
|======================================================================| 100%forecast_h$pred_automl <- h2o.predict(autoML1@leader, forecast_h)##
|
| | 0%
|
|======================================================================| 100%final_forecast <- as.data.frame(forecast_h)
plot_ly(x = df$date, y = df$y,
type = "scatter",
mode = "line",
name = "Actual") %>%
add_lines(x = final_forecast$date, y = final_forecast$pred_rf, name = "Random Forest") %>%
add_lines(x = final_forecast$date, y = final_forecast$pred_gbm, name = "GBM") %>%
add_lines(x = final_forecast$date, y = final_forecast$pred_automl, name = "Auto ML") %>%
layout(title = "Total Vehicle Sales - Final Forecast",
yaxis = list(title = "Thousands of Units", range = c(1100, 1750)),
xaxis = list(title = "Month", range = c(as.Date("2016-01-01"), as.Date("2020-01-01"))))