27
$\begingroup$

I am trying to do time series analysis and am new to this field. I have daily count of an event from 2006-2009 and I want to fit a time series model to it. Here is the progress that I have made:

timeSeriesObj = ts(x,start=c(2006,1,1),frequency=365.25)
plot.ts(timeSeriesObj)

The resulting plot I get is:

Time Series Plot

In order to verify whether there is seasonality and trend in the data or not, I follow the steps mentioned in this post :

ets(x)
fit <- tbats(x)
seasonal <- !is.null(fit$seasonal)
seasonal

and in Rob J Hyndman's blog:

library(fma)
fit1 <- ets(x)
fit2 <- ets(x,model="ANN")

deviance <- 2*c(logLik(fit1) - logLik(fit2))
df <- attributes(logLik(fit1))$df - attributes(logLik(fit2))$df 
#P value
1-pchisq(deviance,df)

Both cases indicate that there is no seasonality.

When I plot the ACF & PACF of the series, here is what I get:

ACF PACF

My questions are:

  1. Is this the way to handle daily time series data? This page suggests that I should be looking at both weekly and annual patterns but the approach is not clear to me.

  2. I do not know how to proceed once I have the ACF and PACF plots.

  3. Can I simply use the auto.arima function?

    fit <- arima(myts, order=c(p, d, q)

*****Updated Auto.Arima results******

When i change the frequency of the data to 7 according to Rob Hyndman's comments here, auto.arima selects a seasonal ARIMA model and outputs:

Series: timeSeriesObj 
ARIMA(1,1,2)(1,0,1)[7]                    

Coefficients:
       ar1      ma1     ma2    sar1     sma1
      0.89  -1.7877  0.7892  0.9870  -0.9278
s.e.   NaN      NaN     NaN  0.0061   0.0162

sigma^2 estimated as 21.72:  log likelihood=-4319.23
AIC=8650.46   AICc=8650.52   BIC=8682.18 

******Updated Seasonality Check******

When I test seasonality with frequency 7, it outputs True but with seasonality 365.25, it outputs false. Is this enough to conclude a lack of yearly seasonality?

timeSeriesObj = ts(x,start=c(2006,1,1),frequency=7)
fit <- tbats(timeSeriesObj)
seasonal <- !is.null(fit$seasonal)
seasonal

returns:

True

while

timeSeriesObj = ts(x,start=c(2006,1,1),frequency=365.25)
fit <- tbats(timeSeriesObj)
seasonal <- !is.null(fit$seasonal)
seasonal

returns:

False
$\endgroup$
8
  • $\begingroup$ What output does str(x) yield? $\endgroup$ Mar 31, 2015 at 11:22
  • $\begingroup$ It yields num [1:1460] 17 12 12 17 13 14 14 5 12 21 ... $\endgroup$ Mar 31, 2015 at 18:15
  • $\begingroup$ Could you post the data? $\endgroup$
    – forecaster
    Mar 31, 2015 at 22:41
  • $\begingroup$ Unfortunately, I can not. $\endgroup$ Apr 1, 2015 at 4:44
  • 3
    $\begingroup$ OK, I would rely on visual inspection and domain knowledge in addition to statistical tests for detecting seasonality. To your question on whether arima can handle multiple seasonality - Sure ARIMA will handle any type of seasonality, R simple doesn't have the capability to handle it. I would look for commercial solutions if there is high inventory/manufacturing cost involved for the product that you are trying to forecast. R has severe limitations for forecasting task like yours. Look at questions on daily forecasting else where in this site. $\endgroup$
    – forecaster
    Apr 1, 2015 at 14:27

3 Answers 3

27
$\begingroup$

Your ACF and PACF indicate that you at least have weekly seasonality, which is shown by the peaks at lags 7, 14, 21 and so forth.

You may also have yearly seasonality, although it's not obvious from your time series.

Your best bet, given potentially multiple seasonalities, may be a tbats model, which explicitly models multiple types of seasonality. Load the forecast package:

library(forecast)

Your output from str(x) indicates that x does not yet carry information about potentially having multiple seasonalities. Look at ?tbats, and compare the output of str(taylor). Assign the seasonalities:

x.msts <- msts(x,seasonal.periods=c(7,365.25))

Now you can fit a tbats model. (Be patient, this may take a while.)

model <- tbats(x.msts)

Finally, you can forecast and plot:

plot(forecast(model,h=100))

You should not use arima() or auto.arima(), since these can only handle a single type of seasonality: either weekly or yearly. Don't ask me what auto.arima() would do on your data. It may pick one of the seasonalities, or it may disregard them altogether.


EDIT to answer additional questions from a comment:

  1. How can I check whether the data has a yearly seasonality or not? Can I create another series of total number of events per month and use its ACF to decide this?

Calculating a model on monthly data might be a possibility. Then you could, e.g., compare AICs between models with and without seasonality.

However, I'd rather use a holdout sample to assess forecasting models. Hold out the last 100 data points. Fit a model with yearly and weekly seasonality to the rest of the data (like above), then fit one with only weekly seasonality, e.g., using auto.arima() on a ts with frequency=7. Forecast using both models into the holdout period. Check which one has a lower error, using MAE, MSE or whatever is most relevant to your loss function. If there is little difference between errors, go with the simpler model; otherwise, use the one with the lower error.

The proof of the pudding is in the eating, and the proof of the time series model is in the forecasting.

To improve matters, don't use a single holdout sample (which may be misleading, given the uptick at the end of your series), but use rolling origin forecasts, which is also known as "time series cross-validation". (I very much recommend that entire free online forecasting textbook.

  1. So Seasonal ARIMA models cannot usually handle multiple seasonalities? Is it a property of the model itself or is it just the way the functions in R are written?

Standard ARIMA models handle seasonality by seasonal differencing. For seasonal monthly data, you would not model the raw time series, but the time series of differences between March 2015 and March 2014, between February 2015 and February 2014 and so forth. (To get forecasts on the original scale, you'd of course need to undifference again.)

There is no immediately obvious way to extend this idea to multiple seasonalities.

Of course, you can do something using ARIMAX, e.g., by including monthly dummies to model the yearly seasonality, then model residuals using weekly seasonal ARIMA. If you want to do this in R, use ts(x,frequency=7), create a matrix of monthly dummies and feed that into the xreg parameter of auto.arima().

I don't recall any publication that specifically extends ARIMA to multiple seasonalities, although I'm sure somebody has done something along the lines in my previous paragraph.

$\endgroup$
5
  • $\begingroup$ Here are a few questions that I have based on your answer: 1. How can I check whether the data has a yearly seasonality or not? Can i create another series of total number of events per month and use its acf to decide this? 2. So Seasonal ARIMA models cannot usually handle multiple seasonalities? Is it a property of the model itself or is it just the way the functions in R are written? $\endgroup$ Apr 1, 2015 at 2:38
  • $\begingroup$ Updated post with results from Auto.Arima with weekly seasonality $\endgroup$ Apr 1, 2015 at 2:50
  • 1
    $\begingroup$ @StephanKolassa, I found this article from AT&T a while back that uses multiple seasonal ARIMA. Model such as the one in article is not possible in R, since R does not have the capability to handle multiseasonal ARIMA. $\endgroup$
    – forecaster
    Apr 1, 2015 at 14:18
  • 2
    $\begingroup$ @forecaster: cool, thanks! It seems like they do double differencing in equation 3.1. I'm kind of concerned about losing a lot of data that way. Unfortunately, they don't compare their results to a simple benchmark, e.g., last week's demands. What I like is how they also investigate forecast combinations between their DSARIMA and a tbats-like model. $\endgroup$ Apr 1, 2015 at 14:28
  • 3
    $\begingroup$ @StephanKolassa I agree, I'm a big proponent of using naive forecast as a benchmark as supported in Principles of Forecasting by Armstrong and only add complexity if it improves accuracy. I ended up in top 2 percentile in few kaggle competitions by using naive methods in my ensembles. $\endgroup$
    – forecaster
    Apr 1, 2015 at 14:33
4
$\begingroup$

The best way to decompose seasonal data using existing R packages is ceemdan() in Rlibeemd. This technique extracts seasonality of multiple periods. The defaults work well. It uses the Hilbert-Huang transform instead of the Fourier transform. The Fourier transform has a severe drawback in that it can only handle stationary, linear data when most series of interest are neither. For example, the random walk y_t = y_{t-1} + e_t is the simplest random walk and frequently encountered. Other methods hold the amplitude of seasonal variation fixed when it often varies in practice.

$\endgroup$
2
  • 1
    $\begingroup$ Please see the basic MathJax tutorial & reference for help on putting mathematics into answers. Please also offer some justification of the claim of "best" (or consider modifying the claim) - it would have to be at least as good as every other option, not just most of them. $\endgroup$
    – Glen_b
    May 25, 2017 at 0:52
  • 2
    $\begingroup$ It might be worth mentioning that this package is on CRAN $\endgroup$
    – Glen_b
    May 25, 2017 at 1:00
3
$\begingroup$

The questions you raise have been dealt with in R Time Series Forecasting: Questions regarding my output . Please look carefully at my detailed answer and all the comments in the discussion including those to the original question as I believe they are relevant to your problem. You might actually take the data that was provided in the post and use it as a teaching moment for yourself. Use the entire discussion as a primer for what you should do.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks a lot! I will use the data from that post to try things out. $\endgroup$ Apr 1, 2015 at 16:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.