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](http://stats.stackexchange.com/questions/47729/two-seasonal-periods-in-arima-using-r) 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](http://robjhyndman.com/hyndsight/longseasonality/) 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](http://www.inside-r.org/packages/cran/forecast/docs/fourier). 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](http://robjhyndman.com/hyndsight/forecasting-weekly-data/) 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](http://robjhyndman.com/hyndsight/longseasonality/) 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)
    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 `arimod`this is the result of the Arima fitting :

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

It doesn't even take into consideration the seasonal part of the model, and 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.

**Edit 2 :** Also I can't seem to be able to add another exogenous variable to the Arima function. I need to use `temperature` - probably as a lead - to fit the `SARIMAX` model, but as soon as I write this :

    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), tmp[1:n]))
    plot(forecast(fit, h=14*m, xreg=cbind(fourier(n+1:(14*m),4,m),fourier(n+1:(14*m),4,1008), tmp[n+1:(14*m)])))

Nothing is plotted besides the initial data set. There is no forecast while without `tmp` as an `xreg` I still get some results.