It is suggested to use auto.arima with xreg in regression with ARIMA errors. Especially when dealing with multiseasonality with regressors, we can use Fourier terms plus dummies and complete with ARIMA errors.
However, auto.arima assumes homoskedasticity. After running auto.arima and determining optimal orders of ARIMA and the number of harmonics, I have tried to run the model step by step.
The regression model with Fourier terms and dummies shows heteroskedasticity. Is the model still available? Should I necessarily correct for heteroskedasticity?
clearly shows a major intercept change at or around period 270 and another intercept change at 65. The analysis time,important to the OP was 21 seconds on a two year old dell portable. The final model required Generalized Least Squares as the error variance changed dramatically at 
week 42 of 2011 . The equation in algebraic form is presented here 
and here 
. Partially presented again here 
with all coefficients being statistically significant at alpha = .05 . The actual and cleansed series highlight the unusual activity 
The Actual/Fit and Forecast plot is 
with the list of forecasts here 
and graphically 
....Looks like the flag ..... If you have any particular detailed questions please contact me offline ... as I don't know how to set up a chat room ....
Best Answer
Yes, it is very important to deal with non constant variance. Ruey Tsay published a paper on this called Outliers, level shifts, and variance changes in time series.