Increasing the accuracy of tbats() forecasts by factoring for correlations between different time-series? This question builds on my previous question Forecasting Hourly Time Series based on previous weeks and same period in previous year/s. My project is to forecast the number of ~400 different types of events expected in each hourly interval with enough accuracy for staffing decisions to be made.
Based on my knowledge of the data I know that each interval is related to the same hour band from the previous few weeks and the same time in the previous few years. Thanks to a comment by Rob Hyndman I am now using tbats() to forecast with mixed results. When comparing the forecasted data to the actual data the monthly totals are consistently within 1-3%.
However, when I compare the forecast to the actual for individual intervals the results are not very reliable at all. I have calculated the difference between the actual and the forecast as a percentage of the forecast for each interval and get an interquartile range of 50% to 150% with a mean difference of ~70%. This level of accuracy is unacceptable for what I need to use the data to do.
I am pretty certain that there are correlations between the frequency of different types of events and some measurable environmental factors. Is there an easy way to feed R a time series for the count of each event type as well as some environmental factors and have it find the correlations and create a forecast?
I am not trying to be lazy, the forecast tool is going to be automated and needs to be able to run without human input.
The method I am currently using is:
data <- scan("data.csv")
fcast <- forecast(tbats(msts(data, seasonal.periods=c(168,8766))),1464)

The csv is a single column containing an hourly count of a specific event type over 2 years.
 A: I would try seasonally adjusted naive methods first, which you could easily employ. Naive method means what you did last year is what you are going to do this year. Of course you have to adjust for Holidays and special events such as day of the week. Before venturing into complex methods I would try this first. It is very hard to beat naive method for short term forecasting.
Looks like you are trying to solve a complex forecasting problem.
There are techniques such as Triple seasonal exponential smoothing by Taylor or by Shen and Huang that you could employ. All the aforementioned methods do not have off the shelf implementation in any software, so you might want to code this yourself.
Having said that, here is what I would do:


*

*If I did not know R/Statistics (I'm not an expert in either one of
these) I would hire a consultant/statistician who    is familiar
with forecasting staffing problems. I would justify these to my
project manager or your supervisor by showing him/her the cost
saving the company would do by using this approach and improving
accuracy of the forecast. If this is done for staffing 400 people, I
would like an expert advice for making a critical decision like
this.

*There are very specialized softwares that does forecasting say for example call    center. I would use those specialized softwares
to solve the problem, of course with proper training and I would use them only if the
software beats my naive bench mark.

A: 
However, when I compare the forecast to the actual for individual intervals the results are not very reliable at all. I have calculated the difference between the actual and the forecast as a percentage of the forecast for each interval and get an interquartile range of 50% to 150% with a mean difference of ~70%. This level of accuracy is unacceptable for what I need to use the data to do.

I know I'm late to the party here, but I wanted to add in case others stumble upon this question.
I can't help but wonder if this is due to not taking into account the daily seasonal component as well. In my opinion, the code should be the following:
data <- scan("data.csv")
fcast <- forecast(tbats(msts(data, seasonal.periods=c(24,168,8766))),1464)

It seems to me that not taking into account the hour of the day would yield the possibility of being way off in the individual intervals.
