PA Discussion 4
Cameron Turner
2024-04-02
Predictive Analytics Discussion 4: ARIMA with External Regressors
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 521495 27.9 1155340 61.8 NA 669291 35.8
## Vcells 960590 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: 0x7f96a3ef6588>
## <environment: namespace:grDevices> #make vector containing the packages
packages <- c("ggplot2",
"knitr",
"fpp3",
"forecast",
"dplyr",
"fable")
#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 # Remove unnecessary values from environment
rm(i,packages)
### Data
# Import the data
unem <- read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Discussions/Data/UNRATENSA.csv")
# Adjust variable formats where necessary
unem$DATE <- as.Date(unem$DATE)
#Filter to past 10 years
unem <- unem %>%
filter(DATE >= max(unem$DATE) - years(10))
## Separate into training and testing data
# Training
unem_ts <- subset(unem, DATE <= max(unem$DATE) - years(1))
# testing data
unem_test <- subset(unem, DATE > max(unem$DATE) - years(1))
# Convert to time series
unem_ts <- ts(unem_ts$UNRATENSA,
start = c(2014,02),
frequency = 12
)
unem_test <- ts(unem_test$UNRATENSA,
start = c(2023,03),
frequency = 12
) Data Observation: 10-YR Real Interest (Same from Discussion 3)
# Show summary stats
summary(unem)## DATE UNRATENSA
## Min. :2014-02-01 Min. : 3.100
## 1st Qu.:2016-08-01 1st Qu.: 3.800
## Median :2019-02-01 Median : 4.400
## Mean :2019-01-30 Mean : 4.869
## 3rd Qu.:2021-08-01 3rd Qu.: 5.500
## Max. :2024-02-01 Max. :14.400#Show head of table
head(unem)## DATE UNRATENSA
## 1 2014-02-01 7.0
## 2 2014-03-01 6.8
## 3 2014-04-01 5.9
## 4 2014-05-01 6.1
## 5 2014-06-01 6.3
## 6 2014-07-01 6.5# Show plot of data
ggplot(unem, aes(x = DATE, y = UNRATENSA)) +
geom_line() +
labs(x = "Date",
y = "Unemployment Rate (%)",
title = "Unemployment Rate")
Models
Dynamic Regression Model- Fed Funds Rate
The regressor used for this model was the Federal Funds Rate. While the link between the two does not appear to be particularly strong, changes in interest rates are generally expected to have a positive effect on the unemployment rate
# Import data for independent variables
# Real GDP
rgdp <- read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Discussions/Data/GDPC1.csv")
rgdp$DATE <- as.Date(rgdp$DATE)
# Fed Funds Rate
ffr <- read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Discussions/Data/FEDFUNDS.csv")
ffr$DATE <- as.Date(ffr$DATE)
# Show plot of data
ggplot(ffr, aes(x = DATE, y = FEDFUNDS)) +
geom_line() +
labs(x = "Date",
y = "Fed Funds Rate (%)",
title = "Fed Funds Rate")
## Separate into training and testing data
# Training
ffr_ts <- subset(ffr, DATE <= max(ffr$DATE) - years(1))
# testing data
ffr_test <- subset(ffr, DATE > max(ffr$DATE) - years(1))
# Convert to time series
ffr_ts <- ts(ffr_ts$FEDFUNDS,
start = c(2014,02),
frequency = 12
)
ffr_test <- ts(ffr_test$FEDFUNDS,
start = c(2023,03),
frequency = 12
)
# Create the model
dyn_mod <- auto.arima(unem_ts,
xreg = ffr_ts, # dynamic regression variable
seasonal = TRUE)
# forecast ARIMA model
dyn_forecast <- forecast(dyn_mod,
xreg = ffr_test,
h = 12)
autoplot(dyn_mod)
# Display arima model
autoplot(unem_ts)+
autolayer(dyn_forecast)+
labs(x = "Date",
y = "Unemployment Rate",
title = "Dynamic Regression Forecast") +
autolayer(unem_test,
series = "True Values")
# Summarize the model to ge
print(summary(dyn_mod))## Series: unem_ts
## Regression with ARIMA(1,0,0) errors
##
## Coefficients:
## ar1 intercept xreg
## 0.7682 5.6454 -0.5870
## s.e. 0.0628 0.5034 0.2933
##
## sigma^2 = 1.029: log likelihood = -155.17
## AIC=318.34 AICc=318.73 BIC=329.11
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01450816 1.000564 0.4052491 -2.220471 7.536722 0.2582731
## ACF1
## Training set 0.09375257# Evaluate accuracy
# Make combined df
dyn_data <- dyn_forecast$mean
dyn_data <- data.frame(unem_test,dyn_data)
# Residuals
dyn_data$residuals <- dyn_data$unem_test - dyn_data$dyn_data
# MSE (for )
MSE_dyn <- mean(dyn_data$residuals^2)
print(paste("MSE for predictions: ",MSE_dyn))## [1] "MSE for predictions: 0.971907137317841"Regular ARIMA
# Create Arima model
arima_mod <- auto.arima(unem_ts,
seasonal = TRUE)
# forecast arima model
arima_forecast <- forecast(arima_mod, h = 12)
autoplot(arima_mod)
# Display arima model
autoplot(unem_ts)+
autolayer(arima_forecast)+
labs(x = "Date",
y = "Unemployment Rate",
title = "ARIMA Forecast") +
autolayer(unem_test,
series = "True Values")
# Summarize the model to ge
print(summary(arima_mod))## Series: unem_ts
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.8260 5.0360
## s.e. 0.0529 0.5352
##
## sigma^2 = 1.047: log likelihood = -156.72
## AIC=319.43 AICc=319.66 BIC=327.51
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02274046 1.013652 0.4052223 -2.619346 7.545573 0.258256
## ACF1
## Training set 0.09232674# Evaluate accuracy
# Make combined df
arima_data <- arima_forecast$mean
arima_data <- data.frame(unem_test,arima_data)
# Residuals
arima_data$residuals <- arima_data$unem_test - arima_data$arima_data
# MSE (for )
MSE_arima <- mean(arima_data$residuals^2)
print(paste("MSE for predictions: ",MSE_arima))## [1] "MSE for predictions: 0.986894937310403"ETS Model
# Create ets model
ets_mod <- ets(unem_ts)
# forecast ets model
ets_forecast <- forecast(ets_mod, h = 12)
autoplot(ets_mod)
# Display ETS model
autoplot(unem_ts)+
autolayer(ets_forecast)+
labs(x = "Date",
y = "Unemployment Rate %",
title = "ETS Forecast") +
autolayer(unem_test,
series = "True Values")
# Evaluate accuracy
# Make combined df
ets_data <- ets_forecast$mean
ets_data <- data.frame(unem_test,ets_data)
# Residuals
ets_data$residuals <- ets_data$unem_test - ets_data$ets_data
# MSE (for )
MSE_ets <- mean(ets_data$residuals^2)
print(paste("MSE for predictions: ",MSE_ets))## [1] "MSE for predictions: 1.203832873331" # Summarize the model to ge
print(summary(ets_mod))## ETS(M,N,A)
##
## Call:
## ets(y = unem_ts)
##
## Smoothing parameters:
## alpha = 0.9047
## gamma = 2e-04
##
## Initial states:
## l = 7.0562
## s = -0.3052 -0.6855 -0.6295 -0.3089 -0.026 0.5251
## 0.8367 0.86 0.3454 0.9107 -0.9165 -0.6062
##
## sigma: 0.1833
##
## AIC AICc BIC
## 496.8214 501.9827 537.1916
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02617299 1.08099 0.47284 -1.114367 8.956728 0.3013501
## ACF1
## Training set 0.04654831Conclusion
Each model had different strengths and weaknesses. To start off, ETS was the weakest model. While it appeared to do its job with predicting seasonal trends, it had high positive bias at the start and seemed to miss the upward trend of the true values. It also had the highest MSE of the three models, making it the least accurate. Of the two ARIMA models, it appeared that the non-regressor model may be a better predictor. While the dynamic regression model had a lower MSE, it had a consistent negative bias and predicted a negative trend. Therefore, the non-regressor ARIMA model may be a better fit at the cost of marginally less accuracy in the short run. It simply seemed to match the steady positive trend of the true values.