While working on a big data set made of 10-minutes-points of information - i.e. 144 points per day, 1008 per week and 52560 per year - I encountered a few problem in R. The information concerns electricity load on a source substation during the year.

###Multiple seasonality :### The data set clearly shows multiple seasonalities, which are daily, weekly and yearly. From there I understood that R doesn't handle multiple seasonality within the ARIMA modeling functions. I would really like to work with ARIMA models though, because my previous work is based on ARIMA models and I know approximatively how to translate a model into an equation.

###Long seasonality :### Each of the seasonalities is of high value, with the shortest one being the daily seasonality at 144. Unfortunately from the SARIMA general equation which is
$\phi(B)\Phi(B^s)W_t = \theta(B)\Theta(B^s)Z_t$
I guessed that the maximum lag for a given model SARIMA(p,d,q)(P,D,Q)144 is
$max((p+P*144), (q+Q*144))$

I would really like to try and fit models with values of P and/or Q greater than 1, but R doesn't allow me since the maximum supported lag = 350. To do so I found this link which is really interesting and led to new functions in the forecast package by M. Hyndman, called fourier and fourierf which you can find here. But since I am not a specialist in forecasting nor in statistics, I have some difficulties understanding how I can make this work.

The thing is I thing this whole fourier regressors package could help me a lot. From what I understood I could use it to simulate the long-seasonality of my data set, maybe use it to simulate multiple seasonality, and even more it could allow me to introduce exogenous variables - which are the temperature and (public holiday + sundays).
I also tried doing some regression following this example but I couldn't make it work because :

Error in forecast.Arima(bestfit, xreg = fourierf(gas, K = 12, h = 1008)) : 
Number of regressors does not match fitted model

I really hope somebody can help me get a better understanding of these functions. Thanks.

Edit : So I tried my best with the fourier example given here but couldn't figure out how it handles the fitting. Here is the code (I copy-pasted M. Hyndman one and adapted to my data set - unsuccessfully) :

n <- 50000
m <- 144
y <- read.table("auch.txt", skip=1)
fourier <- function(t,terms,period)
{
  n <- length(t)
  X <- matrix(,nrow=n,ncol=2*terms)
  for(i in 1:terms)
  {
    X[,2*i-1] <- sin(2*pi*i*t/period)
    X[,2*i] <- cos(2*pi*i*t/period)
  }
  colnames(X) <- paste(c("S","C"),rep(1:terms,rep(2,terms)),sep="")
  return(X)
}
 
library(forecast)
arimod = fit <- Arima(y[1:n,1], order=c(2,1,5), seasonal=c(1,2,8), xreg=cbind(fourier(1:n,4,m),fourier(1:n,4,1008)))
plot(forecast(fit, h=14*m, xreg=cbind(fourier(n+1:(14*m),4,m), fourier(n+1:(14*m),4,1008))))

So I wanted to "force" the model to be a SARIMA(2,1,5)(1,2,8)[144] but when I type arimodthis is the result of the Arima fitting :

> arimod  
Series: y[1:n, 1] , 
ARIMA(2,1,5)                  
  
sigma^2 estimated as 696895:  log likelihood=-407290.2  
AIC=814628.3   AICc=814628.3   BIC=814840

I don't know much about the range the AIC values can take, but it seems way too high to be a good fitting model right there. I think it all comes down to my misunderstanding of the use of Fourier terms as regressors, but I can't figure out why.