PA Homework 3
Cameron Turner
2024-04-18
Predictive Analytics HW 3: Neural Nets/Ensemble
Setup
# R Environment
# Clear the workspace
rm(list = ls()) # Clear environment
gc() # Clear unused memory## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 521677 27.9 1155860 61.8 NA 669291 35.8
## Vcells 961919 7.4 8388608 64.0 16384 1840004 14.1 cat("\f") # Clear the console dev.off # Clear the charts## function (which = dev.cur())
## {
## if (which == 1)
## stop("cannot shut down device 1 (the null device)")
## .External(C_devoff, as.integer(which))
## dev.cur()
## }
## <bytecode: 0x7fb513f2e580>
## <environment: namespace:grDevices> #make vector containing the packages
packages <- c("ggplot2",
"knitr",
"fpp3",
"forecast",
"dplyr",
"fable",
"stargazer",
"gridExtra")
#Install packages
for (i in 1:length(packages)) {# Iterate through each package in the 'packages' vector.
if (!packages[i] %in% rownames(installed.packages())) {#check if package is installed
install.packages(packages[i], #install package if it isn't already
repos = "http://cran.rstudio.com/",
dependencies = TRUE)}
library(packages[i], character.only = TRUE)} # Load the package into the R session.## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──## ✔ tibble 3.2.1 ✔ tsibbledata 0.4.1
## ✔ dplyr 1.1.1 ✔ feasts 0.3.2
## ✔ tidyr 1.3.0 ✔ fable 0.3.3
## ✔ lubridate 1.9.2 ✔ fabletools 0.4.1
## ✔ tsibble 1.1.4## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine # Remove unnecessary values from environment
rm(i,packages)
### Data
# Import vehicle sales data
core_df <- read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Homework/Data/VSales_Unem.csv")
# Adjust variable formats where necessary
core_df$DATE <- as.Date(core_df$DATE, format = "%m/%d/%y")
core_df$VSales <- as.numeric(core_df$VSales)
## Separate into training and testing data
# Training
train <- subset(core_df, DATE <= max(core_df$DATE) - years(1))
# Testing data
test <- subset(core_df, DATE > max(core_df$DATE) - years(1))
# Convert to time series
vsales_train <- ts(train$VSales,
start = c(2019,03),
frequency = 12
)
Unem_train <- ts(train$Unem,
start = c(2019,03),
frequency = 12
)
vsales_test <- ts(test$VSales,
start = c(2023,04),
frequency = 12
)
Unem_test <- ts(test$Unem,
start = c(2023,04),
frequency = 12
)Data Observation
Download sources:
# Show plot of data
ggplot(core_df, aes(x = DATE)) +
geom_line(aes(y = VSales)) +
labs(x = "Date",
y = "Total Sales ($)",
title = "Total Motor Vehicle and Parts Dealer Sales")
Modeling
Neural Network Forecast
## Display sales and Unem Comparison
# Plot for Unem
plot_Unem <- ggplot(core_df, aes(x = DATE)) +
geom_line(aes(y = Unem, color = "Unem")) +
labs(x = "Date",
y = "Percent",
title = "Unem",
color = "Series") +
scale_color_manual(values = c("Unem" = "blue"))
# Plot for VSales
plot_vsales <- ggplot(core_df, aes(x = DATE)) +
geom_line(aes(y = VSales, color = "VSales")) +
labs(x = "Date",
y = "Total Sales ($)",
title = "Total Motor Vehicle and Parts Dealer Sales",
color = "Series") +
scale_color_manual(values = c("VSales" = "black"))
# Combine plots into a grid
grid.arrange(plot_Unem, plot_vsales, ncol = 1)
# Clear up R environment
rm(plot_Unem,plot_vsales)
## Develop the model
# Create the model
nnet_mod <- nnetar(vsales_train,
xreg = Unem_train, # neural network variable
seasonal = TRUE)
# forecast Neural Net model
nnet_forecast <- forecast(nnet_mod,
xreg = Unem_test,
h = 12)
# Display arima model
autoplot(vsales_train)+
autolayer(nnet_forecast)+
labs(x = "Date",
y = "Total Sales ($)",
title = "Neural Network Forecast") +
autolayer(vsales_test,
series = "True Values")
# Summarize the model to ge
print(summary(nnet_mod))## Length Class Mode
## x 49 ts numeric
## m 1 -none- numeric
## p 1 -none- numeric
## P 1 -none- numeric
## scalex 2 -none- list
## scalexreg 2 -none- list
## size 1 -none- numeric
## xreg 49 -none- numeric
## subset 49 -none- numeric
## model 20 nnetarmodels list
## nnetargs 1 -none- list
## fitted 49 ts numeric
## residuals 49 ts numeric
## lags 2 -none- numeric
## series 1 -none- character
## method 1 -none- character
## call 4 -none- call# Evaluate accuracy
# Make combined df
nnet_data <- nnet_forecast$mean
nnet_data <- data.frame(vsales_test,nnet_data)
# Residuals
nnet_data$nnet_residuals <- nnet_data$vsales_test - nnet_data$nnet_dataARIMA
# Create arima model
arima_mod <- auto.arima(vsales_train)
# forecast arima model
arima_forecast <- forecast(arima_mod, h = 12)
autoplot(arima_mod)
# Display ETS model
autoplot(vsales_train)+
autolayer(arima_forecast)+
labs(x = "Date",
y = "Total Sales ($)",
title = "ARIMA Forecast") +
autolayer(vsales_test,
series = "True Values")
# Evaluate accuracy
# Make combined df
arima_data <- arima_forecast$mean
arima_data <- data.frame(vsales_test,arima_data)
# Residuals
arima_data$arima_residuals <- arima_data$vsales_test - arima_data$arima_data
# Summarize the model
print(summary(arima_mod))## Series: vsales_train
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.5455
## s.e. 0.1344
##
## sigma^2 = 33508: log likelihood = -317.85
## AIC=639.7 AICc=639.96 BIC=643.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -14.13849 179.2768 141.2798 -2.651007 11.53997 0.5633353
## ACF1
## Training set -0.01461126ETS
# Create ets model
ets_mod <- ets(vsales_train)
# forecast ets model
ets_forecast <- forecast(ets_mod, h = 12)
autoplot(ets_mod)
# Display ETS model
autoplot(vsales_train)+
autolayer(ets_forecast)+
labs(x = "Date",
y = "Total Sales ($)",
title = "ETS Forecast") +
autolayer(vsales_test,
series = "True Values")
# Evaluate accuracy
# Make combined df
ets_data <- ets_forecast$mean
ets_data <- data.frame(vsales_test,ets_data)
# Residuals
ets_data$ets_residuals <- ets_data$vsales_test - ets_data$ets_dataEnsemble
# Create ensemble (combined forecast model)
combined_forecast <- (nnet_forecast$mean + arima_forecast$mean + ets_forecast$mean)/3
ens_data <- combined_forecast
ens_data <- data.frame(vsales_test,ens_data)
# Residuals
ens_data$ens_residuals <- ens_data$vsales_test -ens_data$ens_data
# Summarize the model
print(summary(ets_mod))## ETS(M,N,N)
##
## Call:
## ets(y = vsales_train)
##
## Smoothing parameters:
## alpha = 0.4022
##
## Initial states:
## l = 1529.2437
##
## sigma: 0.1399
##
## AIC AICc BIC
## 704.1132 704.6465 709.7886
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -12.9365 179.4317 142.6544 -2.661199 11.63865 0.5688163 0.02439772 summary(combined_forecast)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1204 1235 1238 1240 1246 1269Model Comparison
### Display graph of forecasts
autoplot(vsales_train,
ylab = "Total Sales ($)",
size = 0.75) +
ggtitle("Sales Forecasts") +
autolayer(vsales_test,
series = "True Values",
size = 0.75) +
autolayer(nnet_forecast$mean,
series = "Neural Network") +
autolayer(arima_forecast$mean,
series = "ARIMA") +
autolayer(ets_forecast$mean,
series = "ETS") +
autolayer(combined_forecast,
series = "Ensemble") +
scale_colour_manual(values = c("red", "blue", "purple", "green", "darkgreen"))
### Show tables comparing model statistics
# training
# nnet
nnet_train <- data.frame(accuracy(nnet_mod))
rownames(nnet_train) <- "Neural Network"
# ARIMA
arima_train <- data.frame(accuracy(arima_mod))
rownames(arima_train) <- "ARIMA"
# ETS
ets_train <- data.frame(accuracy(ets_mod))
rownames(ets_train) <- "ETS"
#Merge
stats_table_train <- rbind(nnet_train,arima_train,ets_train)
kable(stats_table_train,
caption = "**Stats for Training Data**")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| Neural Network | -0.004331 | 116.2405 | 89.29058 | -1.208651 | 7.348981 | 0.3560348 | -0.0232965 |
| ARIMA | -14.138491 | 179.2768 | 141.27985 | -2.651007 | 11.539966 | 0.5633353 | -0.0146113 |
| ETS | -12.936498 | 179.4317 | 142.65442 | -2.661198 | 11.638649 | 0.5688163 | 0.0243977 |
# testing
stats_table <- data.frame(nnet = nnet_data$nnet_residuals,
arima = arima_data$arima_residuals,
ets = ets_data$ets_residuals,
ens = ens_data$ens_residuals)
# Stats for nnet
nnet_table <- data.frame(
RMSE = sqrt(mean(nnet_data$nnet_residuals^2)),
MAE = mean(abs(nnet_data$nnet_residuals)),
MAPE = mean(abs((nnet_data$nnet_residuals / nnet_data$vsales_test) * 100)),
MASE = mean(abs(nnet_data$nnet_residuals) / mean(abs(diff(nnet_data$vsales_test)))),
ME = mean(nnet_data$nnet_residuals),
MPE = mean((nnet_data$nnet_residuals / nnet_data$vsales_test) * 100)
)
rownames(nnet_table) <- "Neural Network"
#Stats for arima
arima_table <- data.frame(
RMSE = sqrt(mean(arima_data$arima_residuals^2)),
MAE = mean(abs(arima_data$arima_residuals)),
MAPE = mean(abs((arima_data$arima_residuals / arima_data$vsales_test) * 100)),
MASE = mean(abs(arima_data$arima_residuals) / mean(abs(diff(arima_data$vsales_test)))),
ME = mean(arima_data$arima_residuals),
MPE = mean((arima_data$arima_residuals / arima_data$vsales_test) * 100)
)
rownames(arima_table) <- "ARIMA"
#Stats for ets
ets_table <- data.frame(
RMSE = sqrt(mean(ets_data$ets_residuals^2)),
MAE = mean(abs(ets_data$ets_residuals)),
MAPE = mean(abs((ets_data$ets_residuals / ets_data$vsales_test) * 100)),
MASE = mean(abs(ets_data$ets_residuals) / mean(abs(diff(ets_data$vsales_test)))),
ME = mean(ets_data$ets_residuals),
MPE = mean((ets_data$ets_residuals / ets_data$vsales_test) * 100)
)
rownames(ets_table) <- "ETS"
#Stats for Ens
ens_table <- data.frame(
RMSE = sqrt(mean(ens_data$ens_residuals^2)),
MAE = mean(abs(ens_data$ens_residuals)),
MAPE = mean(abs((ens_data$ens_residuals / ens_data$vsales_test) * 100)),
MASE = mean(abs(ens_data$ens_residuals) / mean(abs(diff(ens_data$vsales_test)))),
ME = mean(ens_data$ens_residuals),
MPE = mean((ens_data$ens_residuals / ens_data$vsales_test) * 100)
)
rownames(ens_table) <- "Ensemble"
#Merge
stats_table <- rbind(nnet_table,arima_table,ets_table, ens_table)
kable(stats_table,
caption = "**Stats for Testing Data**")| RMSE | MAE | MAPE | MASE | ME | MPE | |
|---|---|---|---|---|---|---|
| Neural Network | 203.0606 | 189.5461 | 13.718068 | 1.5775297 | 189.54608 | 13.718068 |
| ARIMA | 124.2390 | 105.2371 | 7.723513 | 0.8758542 | 63.45111 | 4.059552 |
| ETS | 130.4280 | 109.1916 | 7.962726 | 0.9087657 | 74.84758 | 4.909986 |
| Ensemble | 144.1999 | 127.0296 | 9.160112 | 1.0572258 | 109.28159 | 7.562536 |
Conclusion
The models yielded some interesting results upon examination. For starters, it was expected that the neural network would perform the best and perhaps it would if more variables were included. Some variables to consider may be inflation, real disposable income, and RGDP (although the model would have to be adapted to quarterly data for GDP). While the ARIMA and Ensemble performed better statistically, the neural network managed to somewhat accurately predict the actual trends in total vehicle sales. This success was also reflected in the neural net having the lowest RMSE. However, since model underestimates the true values more than all other models (reflected in the highest MPE), the value of mean absolute error was higher than most other models.