Forecasting data with multiple seasonality

Question

I'm attempting to forecast the number of taxi rides per hour that occur in NYC. I've turned the data into a time series using 24*7 as the frequency:

taxi_ts <- ts(taxi_train, frequency = 24*7)

I then use decompose() to split my ts into seasonal, trend, and random components.

parts<-decompose(train_ts)

• side question if anyone knows, why does this include alot of NA in the data?

I'm assuming now that I would fit a ARIMA model to my trend component to forecast on that.

After I do this though how would I add the seasonal component back on to the predictions?

Answer 1

I don't think that your attempt to use a frequency of 168 will give you the results that are after as it may be too coarse i.e. crude or unrefined .

https://stats.stackexchange.com/search?q=user%3A3382+hourly+data will give you some pointers as to how I think that you you should proceed with hourly data. Essentially daily habits can impact hourly responses/values.

I have have been routinely able to implement a two-pronged approach where hourly forecasts are developed based on good daily forecasts which are developed based upon :

what day of the week it is

what month your are in

what level changes have occurred

what trend changes have occurred

what days of the month exhibit statistically usual effect

arima structure

what week- of- the month you are in

holiday effects before, on and after

long weekend effects

month-end effects

and possibly/probably rainfall and weather conditions.

With good daily forecasts , my approach to your problem is to construct 24 hourly models using daily totals as an exogenous predictor and identify trends, level shifts , memory structure (arima) at the hourly level.

Answer 2

1) I do not know for sure why you get many NAs in there, but most probably your window is too large. Decompose works by sliding window smoothing - it uses some radius $r$ (usually around 2 times smaller than the seasonality unless you changed it) to go through all values of time series, and sums up $r$ values to the left and right of each time series observation. The first and last $r$ values in the time series will be NA, since the window cannot be estimated ($r$ goes beyond the bounds). In your case you should have around $12 \times 7$ values as NA from each side of time series. However, you still have the necessary component in your $figure$ variable.

2) Take you initial time series, and simply subtract the "figure" from it repeatedly - simply move this window along the time series. This will give you a deseasonalized component (trend + cycle + error). However, if you have multiple seasonalities, you simply either repeat the procedure, but a better way is to use SARIMA - seasonal ARIMA. It will also difference your second season on demand by specifying the large $D$.

3) After you get your SARIMA forecast, simply add your seasonal figure (the one that you subtracted) onto it, step-by-step iteratively again. In case you had 2 seasonal removals and ARIMA, use ARIMA forecast, and first add second seasonal figure iteratively, then add the first seasonal figure iteratively.

However you should try modelling seasonality differently - such a large seasonal figure ($24 \times 7$) will be a bad estimate most probably. Try looking into deseasonalization with smaller period up to 20, or use Fourier series, which is very simple. Check on Hyndman's website: https://robjhyndman.com/publications/complex-seasonality