Forecasting seasonal data I am fairly new to time series modelling. Suppose that I have time series of electricity consumed per day of one company. The series is on workdays only and sometimes there is NA value because of a public holiday. 
I would like to forecast future electricity consumption. I am familiar with the ARIMA modelling in R, however to properly find seasonal component in R I need to (using auto.arima) specificy frequency of the data, which I believe is double seasonal - that is, depends on both workday and month of the year. Because I was unable to fit the ARIMA, I found a post about predicting a time series using linear regression. 
Thus my question is, whether the following is statistically correct approach. To deseason my data, I would start with a linear regression $$y_i = day_i + month_i + \epsilon_i,$$
which gives me estimates $y^{seasonal}_i$ of the seasonality of the data. Then I would model the rest of the series, that is $y_i - y^{seasonal}_i$ via ARIMA using Box-Jenkins for series without seasonal component, so basically to model the trend in the data.
The forecasting would be simply to forecast from both of the models and sum it. Is my reasoning mathematically correct? Can it be used for correlated data? 
 A: Strange, because the Box-Jenkins method is precisely designed to find seasonalites. 
If you know the seasonality, you do not need a Fourier Transform. Just make a regression $y \approx f(workday, month)$. As the variable are qualitative, you need first to transform them in factor (or contrast) as here. Then publish your result so that we can guide you.
Your approach is correct, provided of course the seasonality is additive and the process is stationary (i.e., no trend) and the usual requirements about iid normal errors.
I do not think the correlation between workday and month is will significantly influence the result. You my test it by (1) computing the covariance between workday and month on your actual data and (2) by checking that there is not a large changes in the estimated regression coefficients when a predictor variable is added or deleted.
However, as I said, first check there is no trend. First make a fit $y \approx at + b$ and check $a$ is not  significant. And follow Tukey's advise: make a graph of the raw data as function of time.
