Time Series – Understanding the Period/Cycle of Time Series Data

multiple-seasonalitiespythonstatsmodelstime series

I'm trying to understand the meaning of period/cycle length in time series forecasting. Some functions, such as seasonal_decompose and STL (Python statsmodels package) or models like SARIMA have a period or cycle parameter that indicates 'the period of the series' used (period, seasonal, etc). As I understand it, this would be the period after which the seasonal behaviour of the series repeats (according to this answer).

I'm working with web traffic data that was recorded in 1hour intervals over the span of 8-14 months. My layman's fourier analysis seems to indicate that the data's strongest underlying frequencies are: 24h, 335h (14 days), 675h (28 days), 3720h (155 days), 6h and 12h (in order).

Most online examples assume only 1 frequency in the data, usually as a yearly reoccurrence, and there seems to be little mention of special treatment for time series that have more than a yearly seasonal component. Is such special attention to multiple seasonalities commonly practiced? Is there pre-built code specifically designed for this or are 'single' period functions adapted in some way?

Best Answer

First off, your understanding of period/cycle length is correct.

For hourly web traffic, I would expect the strongest seasonalities for cycles of length 24h (daily seasonality) and 168h (weekly seasonality - typically weekdays and weekends differ strongly). For instance, here is traffic to CV.

You have already found our multiple-seasonalities tag. Its tag wiki contains pointers to literature and commonly used algorithms to model this. My impression is that most such specialized models are implemented in R, rather than Python, but I may be wrong about this.

STL and ARIMA can only deal with a single seasonal period. Depending on your data, you may learn a little something from applying these models separately with the two seasonalities. (E.g., a very simple forecast would be to fit separate ARIMA models for the two seasonalities, then average the forecasts.)

Related Solutions

Solved – very high frequency time series analysis (seconds) and Forecasting (Python/R)

Question #1:

The problem is that in the MLE case, both the Python (statsmodels) and R procedures use state-space models to estimate the likelihood. In an SARIMAX class, the state-space grows linearly (or worse) with the number of seasons (because the state-space form incorporates all intermediate lags too - so if you have a lag at 3600, the state-space form also has all the 3599 intermediate lags).

So you now have a couple of issues - first, you're multiplying 3600+ matrices by each other, which is slow. Even worse, state space models need to be initialized and often they are by default initialized using a stationary initialization that requires solving a 3600 linear system. When I tested a 3600 seasonal order, it wasn't even getting past this part.

The R arima function accepts method='CSS' which uses least-squares (conditional MLE instead of full MLE) to solve the problem. Depending on how the arima function works, it could be much better in your case.

In Python, there aren't many good options. The SARIMAX class accepts a conserve_memory option, but if you do that, you can't forecast. To solve the initialization problem, you can call the initialize_approximate_diffuse method to avoid the 3600 linear system solving. However, even in these cases, you'll be multiplying 3600 x 3600 matrices together, which will be quite slow. I would like to update the SARIMAX class to work with sparse matrices (which would solve this problem) but that's probably quite a ways in the future. I don't know of any non-commercial program that implements state space models using sparse matrices.

Question #5:

This was a bug in the statsmodels code. It has been fixed in the repository (see https://github.com/ChadFulton/statsmodels/issues/2)

Solved – R Time series forecasting: Having issues selecting fourier pairs for ARIMA with regressors

You're hitting the wall because you're exhausting limitations of the first fourier transform fourier(1:n,i,m1). As RandomDude correctly pointed out above, # of transforms i should be less than half period (m1).

However, if, with your code, you run 2 cycles -- one for i, and another for j, where j would be # of transforms for the second seasonality cycle fourier(1:n,j,m3), you would still have a lot of room for model improvement.

This is what I've got from your data, even without dummies, only based on AR, MA, and data seasonality:

library(forecast)
y <- msts(ts, c(7,365)) # multiseasonal ts
fit <- auto.arima(y, seasonal=F, xreg=fourier(y, K=c(3,30)))
fit_f <- forecast(fit, xreg= fourierf(y, K=c(3,30), 180), 180)
plot(fit_f)

I suspect the performance will even improve when holidays are added.

Best Answer

Related Solutions

Solved – very high frequency time series analysis (seconds) and Forecasting (Python/R)

Solved – R Time series forecasting: Having issues selecting fourier pairs for ARIMA with regressors

Related Question