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)
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 :
> 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.
