PA Discussion 6: Neural Nets

Setup

### Environment
  # Clear the workspace
  rm(list = ls()) # Clear environment
  gc()            # Clear unused memory
##          used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 521355 27.9    1154940 61.7         NA   669291 35.8
## Vcells 959623  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: 0x7fb4e1ad25d0>
## <environment: namespace:grDevices>
  #make vector containing the packages
  packages <- c("ggplot2",
                "knitr",
                "fpp3",
                "forecast",
                "dplyr",
                "lubridate",
                "grid",
                "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
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
  # Remove unnecessary values from environment
  rm(i,packages)

Data

Source: FRED Economic Data

The core dataset includes quarterly data from the following (as RGDP did not have monthly data:

  • Unemployment Rate, Percent, Not Seasonally Adjusted (UNRATENSA)

  • Personal Consumption Expenditures: Chain-type Price Index, Index 2017=100, Seasonally Adjusted (PCEPI)

  • Real Gross Domestic Product, Billions of Chained 2017 Dollars, Seasonally Adjusted Annual Rate (GDPC1)

### Manipulate Data
  # Download data to use for analysis
  core_df <- 
    read.csv("/Users/cam/Desktop/PREDICTIVE ANALYTICS/Discussions/Data/Neural_Net_Data.csv")
  
  # Filter down to last 10 years
  core_df$DATE <- as.Date(core_df$DATE) # converts to date format
  
  # Rename variables to simpler terms
  core_df <- core_df %>%
    rename(UNEM = UNRATENSA,
           PCE = PCEPI_PC1,
           RGDP = GDPC1_PC1)
  
  # Convert character vars to numeric
  core_df$PCE <- as.numeric(core_df$PCE)
## Warning: NAs introduced by coercion
  core_df$RGDP <- as.numeric(core_df$RGDP)
## Warning: NAs introduced by coercion
  # Drop 2024 quarter 1 data (contains null values)
  core_df <- core_df %>%
    filter(!(year(DATE) == 2024 & quarter(DATE) == 1))
  
  # Filter to past 10 years
  core_df <- core_df %>%
    filter(DATE >= max(core_df$DATE) - years(10))
  
  # Split into training and testing data
  train <- subset(core_df, DATE <= max(core_df$DATE) - years(2))
  test <- subset(core_df, DATE > max(core_df$DATE) - years(2))
  
  ## Convert to time series objects
  unem_train <- ts(train$UNEM,
                   start = c(2013,4),
                   frequency = 4)
  
  unem_test <- ts(test$UNEM,
                  start = c(2022,1),
                  frequency = 4)
  
  PCE_train <- ts(train$PCE,
                   start = c(2013,4),
                   frequency = 4)
  
  PCE_test <- ts(test$PCE,
                  start = c(2022,1),
                  frequency = 4)
  
  RGDP_train <- ts(train$RGDP,
                   start = c(2013,4),
                   frequency = 4)
  
  RGDP_test <- ts(test$RGDP,
                  start = c(2022,1),
                  frequency = 4)

### Display the data
  ggplot() +
    geom_line(data = core_df, aes(x = DATE, y = UNEM, color = "Unemployment Rate"), linewidth = 1) +
    geom_line(data = core_df, aes(x = DATE, y = PCE, color = "PCE")) +
    geom_line(data = core_df, aes(x = DATE, y = RGDP, color = "RGDP Growth Rate")) +
    labs(title = "Quarterly Data (Past 10 Years)",
         x = "Date",
         y = "Percent, Percent Change from 1 Year Ago",
         color = "Variable") +
    scale_color_manual(values = c("darkgreen", "red", "darkblue"))

Model Creation

Neural Net

# Create Neural Net model
  nnet_mod <- nnetar(unem_train,
                     xreg = cbind(PCE_train,RGDP_train))
  
  # forecast ARIMA model
  nnet_forecast <- forecast(nnet_mod,
                            xreg = cbind(PCE_test,RGDP_test),
                            h = 8)
## Warning in forecast.nnetar(nnet_mod, xreg = cbind(PCE_test, RGDP_test), : xreg
## contains different column names from the xreg used in training. Please check
## that the regressors are in the same order.
  # Display stats for model on training data
  summary(nnet_mod)
##           Length Class        Mode     
## x         33     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      66     mts          numeric  
## subset    33     -none-       numeric  
## model     20     nnetarmodels list     
## nnetargs   0     -none-       list     
## fitted    33     ts           numeric  
## residuals 33     ts           numeric  
## lags       2     -none-       numeric  
## series     1     -none-       character
## method     1     -none-       character
## call       3     -none-       call
  # Display Neural Net model
  autoplot(unem_train)+
    autolayer(nnet_forecast)+
    labs(x = "Date",
         y = "Unemployment Rate %",
         title = "Neural Net Forecast") +
    autolayer(unem_test,
              series = "True Values")

# Evaluate accuracy
  # Make combined df
  nnet_data <- nnet_forecast$mean
  nnet_data <- data.frame(unem_test,nnet_data)
  
  # Rename variables to proper terms
  nnet_data <- nnet_data %>%
    rename(True_Vals = unem_test,
           Forecast_Vals = nnet_data)

  # Residuals
  nnet_data$residuals <- nnet_data$True_Vals - nnet_data$Forecast_Vals

ARIMA

# Create ARIMA model
  arima_mod <- auto.arima(unem_train)
  
  # forecast ARIMA model
  arima_forecast <- forecast(arima_mod, h = 8)

  # Display stats for model on training data
  summary(arima_mod)
## Series: unem_train 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.4914  5.2618
## s.e.  0.1499  0.5206
## 
## sigma^2 = 2.607:  log likelihood = -61.74
## AIC=129.48   AICc=130.31   BIC=133.97
## 
## Training set error measures:
##                       ME     RMSE       MAE       MPE     MAPE      MASE
## Training set -0.02704315 1.564911 0.7665598 -5.976724 13.60213 0.4965801
##                    ACF1
## Training set 0.04187761
  # Display ARIMA model
  autoplot(unem_train)+
    autolayer(arima_forecast)+
    labs(x = "Date",
         y = "Unemployment Rate %",
         title = "ARIMA Forecast") +
    autolayer(unem_test,
              series = "True Values")

# Evaluate accuracy
  # Make combined df
  arima_data <- arima_forecast$mean
  arima_data <- data.frame(unem_test,arima_data)
  
  # Rename variables to proper terms
  arima_data <- arima_data %>%
    rename(True_Vals = unem_test,
           Forecast_Vals = arima_data)

  # Residuals
  arima_data$residuals <- arima_data$True_Vals - arima_data$Forecast_Vals

ETS

# Create ets model
  ets_mod <- ets(unem_train)
  
  # forecast ets model
  ets_forecast <- forecast(ets_mod, h = 8)
  
  # Display stats for model on training data
  summary(ets_mod)
## ETS(A,N,N) 
## 
## Call:
##  ets(y = unem_train) 
## 
##   Smoothing parameters:
##     alpha = 0.6072 
## 
##   Initial states:
##     l = 6.6523 
## 
##   sigma:  1.779
## 
##      AIC     AICc      BIC 
## 157.3425 158.1701 161.8320 
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set -0.1020116 1.724288 0.7741032 -5.251423 12.22742 0.5014667
##                    ACF1
## Training set 0.08687058
  # Display ETS model
  autoplot(unem_train)+
    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)
  
  # Rename variables to proper terms
  ets_data <- ets_data %>%
    rename(True_Vals = unem_test,
           Forecast_Vals = ets_data)

  # Residuals
  ets_data$residuals <- ets_data$True_Vals - ets_data$Forecast_Vals

Model Comparison/Conclusion

# Display graph of all forecasts
  autoplot(unem_train, 
           ylab = "Percent, Percent change from a year ago",
           size = 0.75) +
    ggtitle("Inflation (PCE) Forecasts") +
    autolayer(unem_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") +
    scale_colour_manual(values = c("red", "blue", "purple", "darkgreen"))

## Create stats table describing forecast performances

  # Stats for nnet
  nnet_table <- data.frame(
    RMSE = sqrt(mean(nnet_data$residuals^2)),
    MAE = mean(abs(nnet_data$residuals)),
    MAPE = mean(abs((nnet_data$residuals / nnet_data$True_Vals) * 100)),
    MASE = mean(abs(nnet_data$residuals) / mean(abs(diff(nnet_data$True_Vals)))),
    ME = mean(nnet_data$residuals),
    MPE = mean((nnet_data$residuals / nnet_data$True_Vals) * 100)
  )
  
  rownames(nnet_table) <- "Neural Network"
  
  #Stats for ARIMA
  arima_table <- data.frame(
    RMSE = sqrt(mean(arima_data$residuals^2)),
    MAE = mean(abs(arima_data$residuals)),
    MAPE = mean(abs((arima_data$residuals / arima_data$True_Vals) * 100)),
    MASE = mean(abs(arima_data$residuals) / mean(abs(diff(arima_data$True_Vals)))),
    ME = mean(arima_data$residuals),
    MPE = mean((arima_data$residuals / arima_data$True_Vals) * 100)
  )
  
  rownames(arima_table) <- "ARIMA"  
  
  #Stats for ETS
  ets_table <- data.frame(
    RMSE = sqrt(mean(ets_data$residuals^2)),
    MAE = mean(abs(ets_data$residuals)),
    MAPE = mean(abs((ets_data$residuals / ets_data$True_Vals) * 100)),
    MASE = mean(abs(ets_data$residuals) / mean(abs(diff(ets_data$True_Vals)))),
    ME = mean(ets_data$residuals),
    MPE = mean((ets_data$residuals / ets_data$True_Vals) * 100)
  )
  
  rownames(ets_table) <- "ETS"
  
  #Merge
  stats_table <- rbind(nnet_table,arima_table,ets_table)
  kable(stats_table, 
        caption = "**Stats for Testing Data**")
Stats for Testing Data
RMSE MAE MAPE MASE ME MPE
Neural Network 1.0861828 0.8792349 24.61737 2.600554 -0.8792349 -24.61737
ARIMA 1.5144492 1.4642001 40.91382 4.330733 -1.4642001 -40.91382
ETS 0.9921087 0.9664736 27.00043 2.858584 -0.9664736 -27.00043

Looking at the results from each of the models, it is apparent that the best performing model was the ETS (A,N,N) model. The reason for this may exist outside the predictive power of a neural network in that the unemployment rate appears to have stabilized around its natural rate. The neural network appears to pick up on the trend of unemployment spiking and gradually coming down. Thus, the shape of the predicted neural net forecast mirrors said trend, Given that the economic activity during covid was fairly abnormal, it may have thrown off the forecast.