Over the year 2017 the

smooth

A lot of changes happened in 2017, and it is hard to mention all of them, but the major ones are:

• Introduction of

  ves()

  function (Vector Exponential Smoothing);
• Decrease of computation time of

  es()

  ,

  ssarima()

  and other functions;
• More flexibility in

  es()

  function;
• Development of model selection mechanism for exogenous variables (this will be covered in one of the next posts in a couple of weeks);
• Development of intermittent data model for

  es()

  ,

  ssarima()

  ,

  ces()

  and

  ges()

  (will also be covered in one of the posts in 2018);
• Introduction of

  auto.ges()

  function, which allows selecting the most appropriate GES model for the data (once again, a post on that is planned);
• Switch to

  roxygen2

  package for the maintenance purposes, which simplified my life substantially.

Not all of these things sound exciting, but any house is built brick by brick and some of those bricks are pretty fundamental. Now the new year has started, and I don’t have any intentions of stopping the development of

smooth

• Intermittent VES model. This is driven by one of the topics I’m working together with John Boylan and Stephan Kolassa;
• Exogenous variables in VES model (also connected with that topic);
• Model selection mechanism in

  ves()

  models;
• Conditional intervals for VES model (currently only the version with independent intervals is implemented);
• State-space VARMA, which is needed for my research with Victoria Grigorieva and Yana Salihova;

• sim.iss()

  function, which should allow generating the occurrence part of intermittent data. This should be helpful for another research I plan to conduct together with John Boylan and Patricia Ramos;
• Vector SMA – the idea proposed by Mikhail Artukhov, this is yet another research topic;
• …And a lot of other cool things and features.

But before starting the work in those areas, I’ve decided to measure the performance of the most recent

smooth

ets()

auto.arima()

forecast

• ES – the default exponential smoothing, ETS(Z,Z,Z) model produced by

  es()

  function from

  smooth

  ;
• ESa – the same ETS model, but without multiplicative trends. This was needed in order to make a proper comparison with

  ets()

  function, which neglects the multiplicative trends by default. The call for this model is

  es(datasets[[i]],"ZXZ")

  . I report two versions of this thing: one from

  smooth

  v2.3.0 from CRAN and the other from

  smooth

  v2.3.1 available on GitHub. The latter has the advanced branch and bound mechanism for these types of models, while the former uses the brute force in model selection;
• SSARIMA – State-space ARIMA implemented in

  auto.ssarima()

  function from

  smooth

  ;
• CES – Complex Exponential Smoothing via

  auto.ces()

  function from

  smooth

  ;
• GES – Generalised Exponential Smoothing via

  auto.ges()

  function from

  smooth

• ETS – ETS model implemented in

  ets()

  function from

  forecast

  package;
• ETSm – ETS model with multiplicative trends, done via

  ets(dataset[[i]]$x, allow.multiplicative.trend=TRUE)

  . This makes the function comparable with the default

  es()

  ;
• ARIMA – ARIMA model constructed using

  auto.arima()

  fucntion.

I used five error measures, some of which were discussed in the post about APEs, and aggregated them over the dataset using mean and median values. Here’s the list:

1. MASE – Mean Absolute Scaled Error;
2. sMAE – scaled Mean Absolute Error;
3. RelMAE – Relative Mean Absolute Error (the benchmark in this case is a simple Naive);
4. sMSE – scaled Mean Squared Error (MSE divided by the squared average actuals in-sample);
5. sCE – scaled Cumulative Error (just the sum of errors divided by the average actuals in-sample);

The following code was used in the experiment, so feel free to use it for whatever purposes:

library(Mcomp)
library(Tcomp)
library(smooth)
library(forecast)

datasets <- c(M1,M3,tourism)

methodsNumber <- 8;
methodsNames <- c("ES","ESa","SSARIMA","CES","GES","ETS","ETSm","ARIMA");
test <- es(datasets[[125]]);
testResults <- array(NA,c(methodsNumber,length(datasets),length(test$accuracy)),
                     dimnames=list(methodsNames, NULL, names(test$accuracy)));

#### ES() ####
j <- 1;
for(i in 1:length(datasets)){
    test <- es(datasets[[i]]);
    testResults[j,i,] <- test$accuracy;
    testTimesFull[i,j] <- test$timeElapsed;
}

#### ES() Additive trends ####
j <- 2;
for(i in 1:length(datasets)){
    test <- es(datasets[[i]], "ZXZ");
    testResults[j,i,] <- test$accuracy;
    testTimesFull[i,j] <- test$timeElapsed;
}

#### SSARIMA() ####
j <- 3;
for(i in 1:length(datasets)){
    test <- auto.ssarima(datasets[[i]]);
    testResults[j,i,] <- test$accuracy;
    testTimesFull[i,j] <- test$timeElapsed;
}

#### CES() ####
j <- 4;
for(i in 1:length(datasets)){
    test <- auto.ces(datasets[[i]]);
    testResults[j,i,] <- test$accuracy;
    testTimesFull[i,j] <- test$timeElapsed;
}

#### GES() ####
j <- 5;
for(i in 1:length(datasets)){
    test <- auto.ges(datasets[[i]]);
    testResults[j,i,] <- test$accuracy;
    testTimesFull[i,j] <- test$timeElapsed;
}

#### ETS() ####
j <- 6;
for(i in 1:length(datasets)){
    testStartTime <- Sys.time();
    test <- ets(datasets[[i]]$x);
    testResults[j,i,] <- Accuracy(datasets[[i]]$xx,forecast(test,h=datasets[[i]]$h)$mean,datasets[[i]]$x);
    testTimesFull[i,j] <- Sys.time() - testStartTime;
}

#### ETS() with multiplicative trends ####
j <- 7;
for(i in 1:length(datasets)){
    testStartTime <- Sys.time();
    test <- ets(datasets[[i]]$x, allow.multiplicative.trend=TRUE);
    testResults[j,i,] <- Accuracy(datasets[[i]]$xx,forecast(test,h=datasets[[i]]$h)$mean,datasets[[i]]$x);
    testTimesFull[i,j] <- Sys.time() - testStartTime;
}

#### ARIMA() ####
j <- 8;
for(i in 1:length(datasets)){
    testStartTime <- Sys.time();
    test <- auto.arima(datasets[[i]]$x);
    testResults[j,i,] <- Accuracy(datasets[[i]]$xx,forecast(test,h=datasets[[i]]$h)$mean,datasets[[i]]$x);
    testTimesFull[i,j] <- Sys.time() - testStartTime;
}

selectedMeasures <- c("MASE","sMAE","RelMAE","sMSE","sCE")

### Overall accuracy
competitionResults <- cbind(round(apply(testResults[,,selectedMeasures],c(1,3),mean),3),
      round(apply(testResults[,,selectedMeasures],c(1,3),median),3))*100

# RelMAE should be aggregated using geometric mean rather than arithmetic
competitionResults[,3] <- round(exp(apply(log(testResults[,,"RelMAE"]),1,mean)),3)*100

competitionResults

### Overall time spent on calculations
colSums(testTimesFull)

All the

smooth

forecast

Accuracy()

smooth

I ended up having 5315 forecasts and I calculated those error measures for the whole dataset and for each of the following categories: yearly, quarterly, monthly and other data. You can find 6 tables with the results at the end of this page (the best value is in green bold). Things to notice are:

1. ES performs similar to ETS, when they go through the similar models. e.g. ESa sometimes outperforms ETS and vice versa, but the difference is usually not significant (see, for example, Table 1. But this can also be observed in the other tables);
2. ETS with multiplicative trends performs worse than the one without them. This holds for both

   es()

   and

   ets()

   ;
3. SSARIMA performs (in terms of accuracy) similar to ARIMA on monthly and quarterly data, but fails on yearly data;
4. CES and GES sometimes outperform other models, but this is not consistent;
5. Branch and bound mechanism of

   es()

   is efficient and improves the performance of the function several times (29 minutes for the calculation of ESa v2.3.1 on the whole dataset vs 91 minutes of ESa v2.3.0) – see Table 6 for the details;
6. The three fastest functions in the comparison are:

   ets()

   with the default parameters values,

   es()

   without multiplicative trends and

   ces()

Finally, I’ve conducted Nemenyi test using

nemenyi()

Overall accuracy of models

The graph above tells us that ETS, ES, ESa, ARIMA and ETSm overall perform similar to each other – they are in one group of models. The other models also perform very close to these, but probably slightly less accurate. In order to see if the difference between the models is indeed significant, we would need to conduct the similar experiment on a larger dataset (M4 maybe?).

Computation time ranking

This graph tells us that

ets()

Concluding this New Year experiment, I would say that

smooth

Happy New Year!