PA Homework 2
Cameron Turner
2024-04-14
Vehicle Sales Forecasting
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 519859 27.8 1150665 61.5 NA 669291 35.8
## Vcells 957445 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: 0x7f7afd2db470>
## <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.csv")
# Adjust variable formats where necessary
core_df$DATE <- as.Date(core_df$DATE, format = "%m/%d/%y")
## 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,01),
frequency = 12
)
cpi_train <- ts(train$CPI,
start = c(2019,01),
frequency = 12
)
vsales_test <- ts(test$VSales,
start = c(2023,03),
frequency = 12
)
cpi_test <- ts(test$CPI,
start = c(2023,03),
frequency = 12
)Data Observation
Download source: https://fred.stlouisfed.org/series/TOTALNSA
# 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
Dynamic Regression Forecast
## Display sales and CPI Comparison
# Plot for CPI
plot_cpi <- ggplot(core_df, aes(x = DATE)) +
geom_line(aes(y = CPI, color = "CPI")) +
labs(x = "Date",
y = "Percent Change from 1 Year Ago (%)",
title = "CPI",
color = "Series") +
scale_color_manual(values = c("CPI" = "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_cpi, plot_vsales, ncol = 1)
## Develop the model
# Create the model
dyn_mod <- auto.arima(vsales_train,
xreg = cpi_train, # dynamic regression variable
seasonal = TRUE)
# forecast Dynamic Reg model
dyn_forecast <- forecast(dyn_mod,
xreg = cpi_test,
h = 12)
autoplot(dyn_mod)
# Display arima model
autoplot(vsales_train)+
autolayer(dyn_forecast)+
labs(x = "Date",
y = "Total Sales ($)",
title = "Dynamic Regression Forecast") +
autolayer(vsales_test,
series = "True Values")
# Summarize the model to ge
print(summary(dyn_mod))## Series: vsales_train
## Regression with ARIMA(0,1,0) errors
##
## Coefficients:
## xreg
## 4190.553
## s.e. 4732.856
##
## sigma^2 = 58329105: log likelihood = -507.12
## AIC=1018.24 AICc=1018.51 BIC=1022.03
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 310.8752 7483.044 4089.279 0.0008513905 3.915328 0.3392926
## ACF1
## Training set 0.01248096# Evaluate accuracy
# Make combined df
dyn_data <- dyn_forecast$mean
dyn_data <- data.frame(vsales_test,dyn_data)
# Residuals
dyn_data$dyn_residuals <- dyn_data$vsales_test - dyn_data$dyn_data
# Coefficient for regression
print(paste("Coefficient for CPI: ",dyn_mod$coef))## [1] "Coefficient for CPI: 4190.55330945023"ARIMA
# 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,0)
##
## sigma^2 = 58052933: log likelihood = -507.51
## AIC=1017.02 AICc=1017.11 BIC=1018.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 658.8856 7542.67 4124.286 0.2361365 3.965505 0.3421971 0.02274084ETS
# 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 <- (dyn_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(A,N,N)
##
## Call:
## ets(y = vsales_train)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 99277.159
##
## sigma: 7698.217
##
## AIC AICc BIC
## 1094.434 1094.956 1100.171
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 657.052 7542.682 4122.447 0.2342388 3.963678 0.3420446 0.02278405 summary(combined_forecast)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 129162 129531 130088 130332 130951 131939Model 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(dyn_forecast$mean,
series = "ARIMA with Regression") +
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
# dyn
dyn_train <- data.frame(accuracy(dyn_mod))
rownames(dyn_train) <- "ARIMA with Regression"
# 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(dyn_train,arima_train,ets_train)
kable(stats_table_train,
caption = "**Stats for Training Data**")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | |
|---|---|---|---|---|---|---|---|
| ARIMA with Regression | 310.8752 | 7483.044 | 4089.279 | 0.0008514 | 3.915329 | 0.3392926 | 0.0124810 |
| ARIMA | 658.8856 | 7542.670 | 4124.286 | 0.2361365 | 3.965505 | 0.3421971 | 0.0227408 |
| ETS | 657.0520 | 7542.682 | 4122.447 | 0.2342388 | 3.963678 | 0.3420446 | 0.0227841 |
# testing
stats_table <- data.frame(dyn = dyn_data$dyn_residuals,
arima = arima_data$arima_residuals,
ets = ets_data$ets_residuals,
ens = ens_data$ens_residuals)
# Stats for dyn
dyn_table <- data.frame(
RMSE = sqrt(mean(dyn_data$dyn_residuals^2)),
MAE = mean(abs(dyn_data$dyn_residuals)),
MAPE = mean(abs((dyn_data$dyn_residuals / dyn_data$vsales_test) * 100)),
MASE = mean(abs(dyn_data$dyn_residuals) / mean(abs(diff(dyn_data$vsales_test)))),
ME = mean(dyn_data$dyn_residuals),
MPE = mean((dyn_data$dyn_residuals / dyn_data$vsales_test) * 100)
)
rownames(dyn_table) <- "Dynamic Regression"
#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(dyn_table,arima_table,ets_table, ens_table)
kable(stats_table,
caption = "**Stats for Testing Data**")| RMSE | MAE | MAPE | MASE | ME | MPE | |
|---|---|---|---|---|---|---|
| Dynamic Regression | 7581.417 | 6793.067 | 5.077087 | 5.676801 | 6431.660 | 4.7986461 |
| ARIMA | 1887.131 | 1735.583 | 1.301170 | 1.450385 | 1049.417 | 0.7740715 |
| ETS | 1886.978 | 1735.446 | 1.301069 | 1.450270 | 1049.142 | 0.7738655 |
| Ensemble | 3663.341 | 3388.684 | 2.534985 | 2.831841 | 2843.406 | 2.1155277 |
Conclusion
After observing the results, it is clear that the worst performing model was the dynamic regression while the ETS (ANN) and ARIMA (0,1,0) more or less performed the same. The likely reason for this is that, while the dynamic regression still uses ARIMA, it is still influenced by an external variable. Since it was determined that sales were positively affected by CPI in the regression model, the decline in CPI during the testing period pulled the forecast for sales down when it more-or-less stayed the same. This tendency to underestimate values can be seen clearly on the plots and is also reflected with the highest mean percent error (MPE) in the testing stats table. This also affects the model’s accuracy as can be seen with root mean squared error (RMSE) and mean absolute error (MAE). As for the ETS and ARIMA models, ETS performs marginally better than ARIMA. The reason for this is that all ETS predicted values are higher than ARIMA by exactly .3. Thus, it underestimates sales a little less than ARIMA and is slightly more accurate. Lastly, the ensemble model was pulled downwards by the dynamic regression and caused it to be less efficient than ETS and ARIMA.