Ways to test accuracy of auto.arima in R

Question

I have around 10000 time series showing one particular metric over 5 hours.

I used auto.arima function

In my previous question, people suggested that I have to use auto.arima for each time series, hold off some of data points and test the prediction with my hold off points.

I am holding off 20% of data points (if you see sample out of 40 I will hold off 8) and then let auto.arima predict. Then I can compare generated 8 values with actual 8 values. But is there a formal way to test accuracy in ARIMA model? Is my approach correcT?

Is there a prebuilt function to test the accuracy of Arima.

Here is the code, I can always use x as my time series.

y=auto.arima(x)
plot(forecast(y,h=30))

Sample time series 1

0.0003748,0.0003929,0.0003653,0.0003557,0.0004463,0.000349,0.0003099,0.0003395,0.0003157,0.0002871,0.0002604,0.0002422,0.0001917,0.0002117,0.0002689

time series 2

0.0003977,0.0003481,0.0002413,0.0002069,0.0002127,0.0002108,0.0002003,0.0002174,0.0002098,0.0002069,0.0001955,0.0001926,0.0002108,0.0002146,0.0002079

Both have 40 points. I can hold off 20% of them (8) and compare after auto.arima predicts. But is there a simpler way I can test accuracy?

Answer 1

Let's work with your two series. They don't have 40 observations, only 15 each, so we will set the holdout length to 4. In addition, I'll pretend they are quarterly (frequency=4) so we can see some "time-seriesey" behavior.

Start by defining a list of time series:

> require(forecast)
> series <- list()
> series[[1]] <- ts(c(0.0003748,0.0003929,0.0003653,0.0003557,0.0004463,0.000349,
    0.0003099,0.0003395,0.0003157,0.0002871,0.0002604,0.0002422,
    0.0001917,0.0002117,0.0002689),frequency=4)
> series[[2]] <-    ts(c(0.0003977,0.0003481,0.0002413,0.0002069,0.0002127,0.0002108,
    0.0002003,0.0002174,0.0002098,0.0002069,0.0001955,0.0001926,
    0.0002108,0.0002146,0.0002079),frequency=4)
> holdout <- 4

Next, calculate forecasts. I'll just use lapply() to run over your list of time series and apply an anonymous function that chops off the last 4 data points, fits a model using auto.arima() and forecasts out 4 data points.

> forecasts <- lapply(series,function(foo) {
    subseries <- ts(head(foo,length(foo)-holdout),start=start(foo),frequency=frequency(foo))
    forecast(auto.arima(subseries),h=holdout)
})

> forecasts

[[1]]
     Point Forecast        Lo 80        Hi 80        Lo 95        Hi 95
3 Q4      0.0002604 0.0001984606 0.0003223394 1.656718e-04 0.0003551282
4 Q1      0.0002604 0.0001728044 0.0003479956 1.264341e-04 0.0003943659
4 Q2      0.0002604 0.0001531177 0.0003676823 9.632594e-05 0.0004244741
4 Q3      0.0002604 0.0001365211 0.0003842789 7.094359e-05 0.0004498564

[[2]]
     Point Forecast         Lo 80        Hi 80         Lo 95        Hi 95
3 Q4      0.0001841  1.378356e-04 0.0002303644  0.0001133447 0.0002548553
4 Q1      0.0001727  6.924969e-05 0.0002761503  0.0000144864 0.0003309136
4 Q2      0.0001613 -1.180548e-05 0.0003344055 -0.0001034420 0.0004260420
4 Q3      0.0001499 -1.035005e-04 0.0004033005 -0.0002376426 0.0005374426

Given that you have 10,000 series, lapply() may run a long time. You may want to use a loop instead, with a progress bar.

The result is a list of forecasts, all with the same horizon. Note that each entry has the correct time period information, even if all your series have different starting points and/or frequencies! This is important, because this will be used internally below to match forecasts and actuals.

Finally, to your questions:

Is my approach correct?

Yes it is, and using a holdout sample is exactly the right way to do this!

Is there a prebuilt function to test the accuracy of Arima?

In fact, there is. It's the accuracy() function in the forecastpackage, which will give you more accuracy measures than you could ever want. I'll use mapply(), but you may again want to use a loop.

> result <- mapply(FUN=accuracy,f=forecasts,x=series,SIMPLIFY=FALSE)
> result
 [[1]]
                        ME         RMSE          MAE        MPE     MAPE
Training set -1.036593e-05 4.608252e-05 3.557953e-05  -4.070836 10.20278
Test set     -3.177500e-05 4.328646e-05 3.602500e-05 -15.798730 17.37924
                  MASE       ACF1 Theil's U
Training set 0.5805517 -0.3401003        NA
Test set     0.5878205 -0.1614330  1.072888

[[2]]
                       ME         RMSE          MAE       MPE     MAPE
Training set 3.430307e-06 3.265457e-05 2.294748e-05  1.904284 10.62380
Test set     3.947500e-05 4.395665e-05 3.947500e-05 18.805557 18.80556
                  MASE       ACF1 Theil's U
Training set 0.4167940 -0.2312260        NA
Test set     0.7169824  0.1870022  4.149249

Note that if you use mapply(), you need SIMPLIFY=FALSE, otherwise mapply() will collapse everything into one big matrix, jumbling in-sample and out-of-sample accuracy measures up.

Finally, you can work with sapply() to extract vectors of accuracy measures you are interested in. For instance, you can extract all RMSEs on test sets like this:

> sapply(result,"[","Test set","RMSE")
[1] 4.328646e-05 4.395665e-05