# Timeseries analysis procedure and methods using R

I am working on a small project where we are trying to predict the prices of commodities (Oil, Aluminium, Tin, etc.) for the next 6 months. I have 12 such variables to predict and I have data from Apr, 2008 - May, 2013.

How should I go about prediction? I have done the following:

• Imported data as a Timeseries dataset
• All variable's seasonality tends to vary with Trend, so I am going to multiplicative model.
• I took log of the variable to convert into additive model
• For each variable decomposed the data using STL

I am planning to use Holt Winters exponential smoothing, ARIMA and neural net to forecast. I split the data as training and testing (80, 20). Planning to choose the model with less MAE, MPE, MAPE and MASE.

Am I doing it right?

Also one question I had was, before passing to ARIMA or neural net should I smooth the data? If yes, using what? The data shows both Seasonality and trend.

EDIT:

Attaching the timeseries plot and data

Year  <- c(2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2009, 2009,
2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2010,
2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,
2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011,
2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012,
2012, 2012, 2013, 2013)
Month <- c(4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2)
Coil  <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000,
29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500,
33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500,
33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700,
40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500,
40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300,
39300, 39300, 39300, 39300, 39800)
coil <- data.frame(Year = Year, Month = Month, Coil = Coil)


EDIT 2: One question, can you please tell me if my data has any seasonality or trend? And also please give me some tips on how to identify them.

• If you're looking at trying to forecast groups of commodities, such as various types of metal (steel A, steel B, steel C, etc.), then it might be worth testing for the existence of cointegration. For example, something like this: Do steel prices move together?. This may provide better 6 month (medium/long term) forecasts than univariate methods, but this is, indeed, a difficult game that you're trying to play. ;-) – Graeme Walsh Mar 4 '15 at 1:06
• AS @GraemeWalsh points out, univariate trend extrapolation might not be ideal for this type of data. There are well established methods in the literature for forecasting oil, steel prices that might be worth exploring. – forecaster Mar 4 '15 at 16:00
• Can you post new edits as a separate question ? Since you have already accepted an answer the new questions might not get the attention it needs. From eyeballing the data I can say none of them have trends or seasonal patterns. As noted in my post below, looks like the downward trend before 2009 is a macro economic phenomena like recession ? – forecaster Mar 4 '15 at 21:19
• @forecaster,@GraemeWalsh : Thank you. I am planning to use cointegration method using ADF tests. – Niranjan Sonachalam Mar 4 '15 at 21:19
• You have provided context in your new question and it makes mosre sense now. So drop prior to 2009 was indeed some macro economic phenomena. In that instance please use random walk method with drift or (arima(0,1,0)+drift – forecaster Mar 4 '15 at 22:25

You should use the forecast package, which supports all of these models (and more) and makes fitting them a snap:

library(forecast)
x <- AirPassengers
mod_arima <- auto.arima(x, ic='aicc', stepwise=FALSE)
mod_exponential <- ets(x, ic='aicc', restrict=FALSE)
mod_neural <- nnetar(x, p=12, size=25)
mod_tbats <- tbats(x, ic='aicc', seasonal.periods=12)
par(mfrow=c(4, 1))
plot(forecast(mod_arima, 12), include=36)
plot(forecast(mod_exponential, 12), include=36)
plot(forecast(mod_neural, 12), include=36)
plot(forecast(mod_tbats, 12), include=36)


I would advise against smoothing the data prior to fitting your model. Your model is inherently going to try to smooth the data, so pre-smoothing just complicates things.

Edit based on new data:

It actually looks like arima is one of the worst models you could chose for this training and test set.

I saved your data to a file call coil.csv, loaded it into R, and split it into a training and test set:

library(forecast)
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12) test_x <- window(x, start=c(2012, 3)) x <- window(x, end=c(2012, 2))  Next I fit a bunch of time series models: arima, exponential smoothing, neural network, tbats, bats, seasonal decomposition, and structural time series: models <- list( mod_arima = auto.arima(x, ic='aicc', stepwise=FALSE), mod_exp = ets(x, ic='aicc', restrict=FALSE), mod_neural = nnetar(x, p=12, size=25), mod_tbats = tbats(x, ic='aicc', seasonal.periods=12), mod_bats = bats(x, ic='aicc', seasonal.periods=12), mod_stl = stlm(x, s.window=12, ic='aicc', robust=TRUE, method='ets'), mod_sts = StructTS(x) )  Then I made some forecasts and compared to the test set. I included a naive forecast that always predicts a flat, horizontal line: forecasts <- lapply(models, forecast, 12) forecasts$naive <- naive(x, 12)
par(mfrow=c(4, 2))
for(f in forecasts){
plot(f)
lines(test_x, col='red')
}


As you can see, the arima model gets the trend wrong, but I kind of like the look of the "Basic Structural Model"

Finally, I measured each model's accuracy on the test set:

acc <- lapply(forecasts, function(f){
accuracy(f, test_x)[2,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
ME    RMSE     MAE   MPE MAPE MASE ACF1 Theil's U
mod_sts     283.15  609.04  514.46  0.69 1.27 0.10 0.77      1.65
mod_bats     65.36  706.93  638.31  0.13 1.59 0.12 0.85      1.96
mod_tbats    65.22  706.92  638.32  0.13 1.59 0.12 0.85      1.96
mod_exp      25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
naive        25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
mod_neural   81.14  853.86  754.61  0.18 1.89 0.14 0.14      2.39
mod_arima   766.51  904.06  766.51  1.90 1.90 0.14 0.73      2.48
mod_stl    -208.74 1166.84 1005.81 -0.52 2.50 0.19 0.32      3.02


The metrics used are described in Hyndman, R.J. and Athanasopoulos, G. (2014) "Forecasting: principles and practice", who also happen to be the authors of the forecast package. I highly recommend you read their text: it's available for free online. The structural time series is the best model by several metrics, including MASE, which is the metric I tend to prefer for model selection.

One final question is: did the structural model get lucky on this test set? One way to assess this is looking at training set errors. Training set errors are less reliable than test set errors (because they can be over-fit), but in this case the structural model still comes out on top:

acc <- lapply(forecasts, function(f){
accuracy(f, test_x)[1,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
ME    RMSE     MAE   MPE MAPE MASE  ACF1 Theil's U
mod_sts      -0.03    0.99    0.71  0.00 0.00 0.00  0.08        NA
mod_neural    3.00 1145.91  839.15 -0.09 2.25 0.16  0.00        NA
mod_exp     -82.74 1915.75 1359.87 -0.33 3.68 0.25  0.06        NA
naive       -86.96 1936.38 1386.96 -0.34 3.75 0.26  0.06        NA
mod_arima  -180.32 1889.56 1393.94 -0.74 3.79 0.26  0.09        NA
mod_stl     -38.12 2158.25 1471.63 -0.22 4.00 0.28 -0.09        NA
mod_bats     57.07 2184.16 1525.28  0.00 4.07 0.29 -0.03        NA
mod_tbats    62.30 2203.54 1531.48  0.01 4.08 0.29 -0.03        NA


(Note that the neural network overfit, performing excellent on the training set and poorly on the test set)

Finally, it would be a good idea to cross-validate all of these models, perhaps by training on 2008-2009/testing on 2010, training on 2008-2010/testing on 2011, training on 2008-2011/testing on 2012, training on 2008-2012/testing on 2013, and averaging errors across all of these time periods. If you wish to go down that route, I have a partially complete package for cross-validating time series models on github that I'd love you to try out and give me feedback/pull requests on:

devtools::install_github('zachmayer/cv.ts')
library(cv.ts)


Edit 2: Lets see if I remember how to use my own package!

First of all, install and load the package from github (see above). Then cross-validate some models (using the full dataset):

library(cv.ts)
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12) ctrl <- tseriesControl(stepSize=1, maxHorizon=12, minObs=36, fixedWindow=TRUE) models <- list() models$arima = cv.ts(
x, auto.arimaForecast, tsControl=ctrl,
ic='aicc', stepwise=FALSE)

models$exp = cv.ts( x, etsForecast, tsControl=ctrl, ic='aicc', restrict=FALSE) models$neural = cv.ts(
x, nnetarForecast, tsControl=ctrl,
nn_p=6, size=5)

models$tbats = cv.ts( x, tbatsForecast, tsControl=ctrl, seasonal.periods=12) models$bats = cv.ts(
x, batsForecast, tsControl=ctrl,
seasonal.periods=12)

models$stl = cv.ts( x, stl.Forecast, tsControl=ctrl, s.window=12, ic='aicc', robust=TRUE, method='ets') models$sts = cv.ts(x, stsForecast, tsControl=ctrl)

models$naive = cv.ts(x, naiveForecast, tsControl=ctrl) models$theta = cv.ts(x, thetaForecast, tsControl=ctrl)


(Note that I reduced the flexibility of the neural network model, to try to help prevent it from overfitting)

Once we've fit the models, we can compare them by MAPE (cv.ts doesn't yet support MASE):

res_overall <- lapply(models, function(x) x$results[13,-1]) res_overall <- Reduce(rbind, res_overall) row.names(res_overall) <- names(models) res_overall <- res_overall[order(res_overall[,'MAPE']),] round(res_overall, 2) ME RMSE MAE MPE MAPE naive 91.40 1126.83 961.18 0.19 2.40 ets 91.56 1127.09 961.35 0.19 2.40 stl -114.59 1661.73 1332.73 -0.29 3.36 neural 5.26 1979.83 1521.83 0.00 3.83 bats 294.01 2087.99 1725.14 0.70 4.32 sts -698.90 3680.71 1901.78 -1.81 4.77 arima -1687.27 2750.49 2199.53 -4.23 5.53 tbats -476.67 2761.44 2428.34 -1.23 6.10  Ouch. It would appear that our structural forecast got lucky. Over the long term, the naive forecast makes the best forecasts, averaged across a 12-month horizon (the arima model is still one of the worst models). Let's compare the models at each of the 12 forecast horizons, and see if any of them ever beat the naive model: library(reshape2) library(ggplot2) res <- lapply(models, function(x) x$results$MAPE[1:12]) res <- data.frame(do.call(cbind, res)) res$horizon <- 1:nrow(res)
res <- melt(res, id.var='horizon', variable.name='model', value.name='MAPE')
res$model <- factor(res$model, levels=row.names(res_overall))
ggplot(res, aes(x=horizon, y=MAPE, col=model)) +
geom_line(size=2) + theme_bw() +
theme(legend.position="top") +
scale_color_manual(values=c(
"#1f78b4", "#ff7f00", "#33a02c", "#6a3d9a",
"#e31a1c", "#b15928", "#a6cee3", "#fdbf6f",
"#b2df8a")
)


Tellingly, the exponential smoothing model is always picking the naive model (the orange line and blue line overlap 100%). In other words, the naive forecast of "next month's coil prices will be the same as this month's coil prices" is more accurate (at almost every forecast horizon) than 7 extremely sophisticated time series models. Unless you have some secret information the coil market doesn't already know, beating the naive coil price forecast is going to be extremely difficult.

It's never the answer anyone wants to hear, but if forecast accuracy is your goal, you should use the most accurate model. Use the naive model.

• It's interesting to look at the differences bt these models. NNAR in particular looks different. Given that this is a famous dataset (& historically old, I believe), is it known which is right & whether one model type outperforms? (Nb, I know relatively little about TS.) – gung Mar 3 '15 at 20:23
• @gung The best way to do this would be to split off a holdout set and test the model. Note that the model that makes the best short-term forecasts may not be the model that makes the best long-term forecasts.... – Zach Mar 3 '15 at 21:00
• Thanks a lot, but I am not getting good forecasts for the above dataset (I think I am missing some important step here). Can you please let me know if I am missing something – Niranjan Sonachalam Mar 4 '15 at 0:28
• @Niranjan Can you tell us/show how you conclude that you are not getting good forecast? – forecaster Mar 4 '15 at 0:32
• @forecaster: Please check here pbrd.co/1DRPRsq . I am new to forecasting. Let me know if you need any specific info. I tried with Arima model. – Niranjan Sonachalam Mar 4 '15 at 0:41

The approach that you have taken is reasonable. If you are new to forecasting, then I would recommend following books:

1. Forecasting methods and applications by Makridakis, Wheelright and Hyndman
2. Forecasting: Principles and practice by Hyndman and Athana­sopou­los.

The first book is a classic which I strongly recommend. The second book is an open source book which you can refer for forecasting methods and how it is applied using R open source software package forecast. Both the books provide good background on the methods that I have used. If you are serious about forecasting, then I would recommend Principles of Forecasting by Armstrong which is collection of tremendous amount of research in forecasting that a practitioners might find it very helpful.

Coming to your specific example on coil, it reminds me of a concept of forecastability which most textbooks often ignore. Some series such as your series simply cannot be forecasted because it is pattern less as it doesn't exhibit trend or seasonal patters or any systematic variation. In that case I would categorize a series as less forecastable. Before venturing into extrapolation methods, I would look at the data and ask the question, is my series forecastable?In this specific example, a simple extrapolation such as random walk forecast which uses the last value of the forecast has been found to be most accurate.

Also one additional comment about neural network: Neural networks are notoriously know to fail in empirical competitions. I would try traditional statitical methods for time series before attempting to use neural networks for time series forecasting tasks.

I attempted to model your data in R's forecast package, hopefully the comments are self explanatory.

coil <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000,
29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500,
33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500,
33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700,
40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500,
40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300,
39300, 39300, 39300, 39300, 39800)

coilts <- ts(coil,start=c(2008,4),frequency=12)

library("forecast")

# Data for modeling
coilts.mod <- window(coilts,end = c(2012,3))

#Data for testing
coil.test <- window(coilts,start=c(2012,4))

# Model using multiple methods - arima, expo smooth, theta, random walk, structural time series

#arima
coil.arima <- forecast(auto.arima(coilts.mod),h=11)

#exponential smoothing
coil.ets <- forecast(ets(coilts.mod),h=11)

#theta
coil.tht <- thetaf(coilts.mod, h=11)

#random walk
coil.rwf <- rwf(coilts.mod, h=11)

#structts
coil.struc <- forecast(StructTS(coilts.mod),h=11)

##accuracy

arm.acc <- accuracy(coil.arima,coil.test)
ets.acc <- accuracy(coil.ets,coil.test)
tht.acc <- accuracy(coil.tht,coil.test)
rwf.acc <- accuracy(coil.rwf,coil.test)
str.acc <- accuracy(coil.struc,coil.test)


Using MAE on the hold out data, I would choose ARIMA for short term forecast (1 - 12 months). for long term, I would rely on random walk forecast. Please note that ARIMA picked a random walk model with drift (0,1,0)+drift which tends to be much more accurate than pure random walk model in these type of problems specifically short term. See below chart. This is based on the accuracy function as shown in the above code.

Specific answers to your specific questions: Also one question I had was, before passing to ARIMA or neural net should I smooth the data? If yes, using what?

• No, Forecasting methods naturally smooths your data to fit model.

The data shows both Seasonality and trend.

• The above data doesnt show trend or seasonality. If you determine that the data exhibits seasonality and trend, then choose an appropriate method.

Practical Tips to improve accuracy:

Combine variety of forecasting methods: - You could try using non extrapolation methods such as forecasting by analogy, judgmental forecasting or other techniques and combine them with your statitical methods to provide accurate predictions. See this article for benefits of combining. I tried combining the above 5 methods, but the prediction were not accurate as individual methods, one possible reason is that individual forecast are similar. You reap the benefits of combining forecast when you combine diverse methods such as statistical and judgmental forecasts.

Detect and Understand Outliers: - Real world data is filled with outliers. Identify and appropriately treat outliers in time series. Recommend reading this post. In looking at your coil data, is the drop prior to 2009 is an outlier ??

Edit

The data appears to be following some type of macro economic trends. My guess is the downward trend seen before 2009 follows slump in economy seen between 2008 - 2009 and start to pick up post 2009. If this is the case, then I would all together avoid using any extrapolation methods and instead rely on sound theory on how these economic trends behave such as the one referenced by @GraemeWalsh.

Hope this helps