# Forecasting data with multiple seasonality

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?

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 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.

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