7.1 Using ndiff function to identify the needed number of differences

ndiffs(employment_training, alpha = 0.05, test = "adf")

## [1] 1

The function calculate the need of just one difference.

7.2 using auto.arima

employment_arima = auto.arima(employment_training)
summary(employment_arima)

## Series: employment_training 
## ARIMA(0,1,0) 
## 
## sigma^2 estimated as 3573478:  log likelihood=-1156
## AIC=2315   AICc=2315   BIC=2317
## 
## Training set error measures:
##                ME RMSE MAE   MPE  MAPE  MASE   ACF1
## Training set 97.1 1883 441 0.062 0.323 0.152 0.0039

The auto.arima() results (ARIMA(0,1,0)) confirmed the initial analysis with ndiffs().

tsdisplay(diff(employment_training))

Applying one difference, the results of the ACF and the PACF are much better, and the graph is showing a stationary data.

The graph also shows a sharp drop in the number of employments around March 2020, probably caused by the COVID-19 pandemic’s lockdowns. However, this is not affecting the stationarity of the chart.

8.1 Forecasting with the arima model and plotting the forecast

employment_arima_fc = forecast(employment_arima, h = length(employment_testing))
autoplot(employment_arima_fc, col = "blue") +
  autolayer(employment_arima_fc$mean, alpha = 0.5, color = "red") 

The mean of the arima model looks like a naive model, always repeating the value of the last period. This is due to having AR(p) and MA(q) igual to zero.

8.2 Calculating the accuracy of the predictions

accuracy(employment_arima_fc$mean, employment_testing)

##            ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## Test set 2758 3287 2758 1.89 1.89 0.767      6.27

The RMSE (3287) of the testing time-series is much higher than the RMSE (1909) of the training time-series, which is something expected because the model was created based on the data of the training time-series.

The forecast data is presenting a linear behavior with zero slope because the model is not using past periods values (p) or error values of past periods forecasts (q) to make the predictions.

I will try to get a model with better RMSE by adding 1 q to the Arima function. The function auto.arima not necessarily find the model with the best RMSE, the function calculates a trade off between accuracy and parsimony based on the Akaike’s Information Criterion (AIC) or/and the Bayesian Information Criterion (BIC).

If our focus is to find the best prediction model regardless of its complexity, we should try other combinations of p, d, and q. I can probably find models with better RMSE, however, they will have worse (higher) AIC and BIC, and may present overfitting issues.

arima_brute_force <-
  function(training_ts,
           testing_ts,
           p_vec = c(0:10),
           d_vec = c(0:10),
           q_vec = c(0:10)) {
    best_model = auto.arima(training_ts)
    best_RMSE = accuracy(best_model)[2]
    for (p in p_vec) {
      for (d in d_vec) {
        for (q in q_vec) {
          err <- FALSE
          tryCatch({
            result <- Arima(training_ts,
                            order = c(p, d, q))
          },
          error = function(e) {
            err <<- TRUE
          },
          warning = function(w) {
          })
          if (!err) {
            r = accuracy(result)
            if (!is.na(accuracy(result)[2])) {
              if (accuracy(result)[2] < best_RMSE) {
                best_model = result
                best_RMSE = accuracy(result)[2]
              }
            }
          }
        }
      }
    }
    result <- best_model
  }

best_model = arima_brute_force(employment_training, employment_testing, c(0:10), c(0:10), c(0:10))
summary(best_model)

## Series: training_ts 
## ARIMA(3,0,8) with non-zero mean 
## 
## Coefficients:

## Warning in sqrt(diag(x$var.coef)): NaNs produced

##         ar1    ar2     ar3     ma1     ma2     ma3    ma4    ma5    ma6    ma7
##       0.999  0.995  -0.997  -0.113  -1.304  -0.025  0.178  0.053  0.059  0.085
## s.e.    NaN    NaN   0.009   0.103   0.107   0.144  0.146  0.142  0.144  0.112
##         ma8    mean
##       0.067  138874
## s.e.  0.112    1680
## 
## sigma^2 estimated as 3375787:  log likelihood=-1159
## AIC=2344   AICc=2347   BIC=2381
## 
## Training set error measures:
##                ME RMSE MAE     MPE  MAPE  MASE     ACF1
## Training set 3.07 1750 489 -0.0177 0.349 0.168 -0.00612

employment_arima_fc_best_model = forecast(best_model, h = length(employment_testing))
autoplot(employment_arima_fc_best_model, col = "blue") +
  autolayer(employment_arima_fc_best_model$mean, alpha = 0.5, color = "red") 

accuracy(employment_arima_fc_best_model$mean, employment_testing)

##            ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## Test set 4548 5364 4548 3.11 3.11 0.767      10.2

10.1 Forecasting with the arimax model and plotting the forecast

employment_arimax_fc = forecast(employment_arimax, xreg = regressors_testing_ts)

## Warning in forecast.forecast_ARIMA(employment_arimax, xreg =
## regressors_testing_ts): xreg contains different column names from the xreg used
## in training. Please check that the regressors are in the same order.

autoplot(employment_arimax_fc, col = "blue") +
  autolayer(employment_arimax_fc$mean, alpha = 0.5, color = "red") 

The graph is also showing some significant improvement.

10.2 Calculating the accuracy of the predictions

accuracy(employment_arimax_fc$mean, employment_testing)

##             ME RMSE  MAE    MPE  MAPE  ACF1 Theil's U
## Test set -1140 1419 1191 -0.788 0.823 0.416      2.74

The RMSE is worse than the RMSE using the training time-series. However, comparing the RMSE of the two testing time-series, the RMSE (1909) of the arimax model is showing a significant improvement in relation to the RMSE (3287) to the arima model

11.1 Forecasting with the NNETS model and plotting the forecast

employment_nnets_fc = forecast(employment_nnets, xreg = regressors_testing_ts)

## Warning in forecast.nnetar(employment_nnets, xreg = regressors_testing_ts): xreg
## contains different column names from the xreg used in training. Please check
## that the regressors are in the same order.

autoplot(employment_ts, col = "blue") +
  autolayer(employment_nnets_fc$mean, color="red") 

11.2 Calculating the accuracy of the predictions

accuracy(employment_nnets_fc, employment_testing)

##                   ME RMSE  MAE      MPE  MAPE   MASE    ACF1 Theil's U
## Training set   -1.03  292  164 -0.00242 0.115 0.0564 -0.0915        NA
## Test set     -864.14 2245 1757 -0.60799 1.215 0.6035  0.7885      4.35

The RMSE (308) of the training set is by far the best one. however, it is presenting a worse RMSE (2512) than the arimax model RMSE (1986). This is a indication of overfitting. If I have to choose one model now I would choose the arimax model with regressors.

12.1 Forecasting with the TBATS model and plotting the forecast

employment_tbats_fc = forecast(employment_tbats, h = length(employment_testing))
autoplot(employment_ts, col = "blue") +
  autolayer(employment_tbats_fc, alpha = 0.5, color = "red") 

12.2 Calculating the accuracy of the predictions

accuracy(employment_tbats_fc$mean, employment_testing)

##            ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## Test set 2751 3282 2751 1.88 1.88 0.767      6.26

The tbats is presenting the worst RMSE for the training (1912) and testing (3276) sets. So I will keep the arimax model with regressors.

13.1 Forecasting with the hybrid model and plotting the forecast

employment_hybrid_fc = forecast(employment_hybrid,
                                h = length(employment_testing),
                                xreg = regressors_testing_ts)

## Warning in forecast.forecast_ARIMA(object$auto.arima, h = h, xreg = xregA, :
## xreg contains different column names from the xreg used in training. Please
## check that the regressors are in the same order.

## Warning in forecast.nnetar(object$nnetar, h = h, xreg = xregN, PI = PI, : xreg
## contains different column names from the xreg used in training. Please check
## that the regressors are in the same order.

autoplot(employment_ts, col = "blue") +
  autolayer(employment_hybrid_fc, alpha = 0.5, color = "red") 

13.2 Calculating the accuracy of the predictions

accuracy(employment_hybrid_fc$mean, employment_testing)

##            ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## Test set 1495 2585 2129 1.01 1.46 0.835      4.93

The hybrid model is the clear winner. The best RMSE (308) for the training set was the nnetar model. However, the nnetar presented low accuracy because of overfitting. The hybrid model had the second-best RMSE (1224) for the training set and the best RMSE (1176) for the testing set, showing no signs of overfitting.