6
$\begingroup$

I've been working on some various time series forecasts and I've begun to notice a trend (pardon the pun) in my analyses. For about 5-7 datasets that I've worked with so far, it would be helpful to allow for multiple seasonal periods along with an option for holiday dummies. I've tried various methods and usually stick with tbats since auto.arima() with regressors has been giving me issues. By this point, it's probably obvious I'm working in R.

Before I get too far, let me give some sample data. Hopefully the following link works: https://gist.github.com/JaredRayWolf/c8cb601dd26ec72a64d0.

This data yields the following time series plot: Time Series Plot The large dips are around Christmas and New Years, however there are also smaller dips around Thanksgiving. In the code below, I name this dataset traindata.

Now, ets and "plain" auto.arima don't look so hot in the long run since they are limited to only one seasonal period (I choose weekly). However for my test set that I held out they performed fairly well for the month's worth of data (with the exception of Labor Day weekend). This being said, forecasting out for a year would be ideal.

I next tried tbats with weekly and yearly seasonal periods. That results in the following forecast: TBATS Forecast

Now this looks pretty good. From the naked eye it looks great at taking into account the weekly and yearly seasonal periods as well as Christmas and New Years effects (since they obviously fall on the same dates each year). It would be best if I could include the holidays (and the days around them) as dummy variables. Hence my attempts at auto.arima with xreg regressors.

For ARIMA with regressors, I've followed Dr. Hyndman's suggestions for the fourier function (given here: http://robjhyndman.com/hyndsight/longseasonality/) as well as his selection of the number of fourier terms (given here: http://robjhyndman.com/hyndsight/forecasting-weekly-data/)

My code is as follows:

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(paste(c("S","C"),rep(1:terms,rep(2,terms)),sep=""),period,sep="_")
  return(X)
}

fcdaysout<-365
m1<-7
m2<-30.4375
m3<-365.25

hol<-cbind(traindata$CPY_HOL, traindata$DAY_BEFORE_CPY_HOL, traindata$DAY_AFTER_CPY_HOL)
hol<-as.matrix(hol)

n <- nrow(traindata)
bestfit <- list(aicc=Inf)
bestk <- 0

for(i in 1:m1)
{
    fake_xreg = cbind(fourier(1:n,i,m1), fourier(1:n,i,m3), hol)
    fit <- auto.arima(traindata$ACTIVE_LOADS, xreg = fake_xreg, max.p=10, max.q=10, max.d=2, stepwise=FALSE, ic="aicc", allowdrift=TRUE)
	if(fit$aicc < bestfit$aicc)
    {
        bestfit <- fit
        bestk <- i
    }
    else
    {
    }
}

k <- bestk
k
##k<-3

xreg<-cbind(fourier(1:n,k,m1), fourier(1:n,k,m3), hol)
xreg<-as.matrix(xreg)

aacov_fit <- auto.arima(traindata$ACTIVE_LOADS, xreg = xreg, max.p=10, max.q=10, max.d=2, stepwise=FALSE, ic="aic", allowdrift=TRUE)
summary(aacov_fit)

Where my issues come in is inside the for loop to determine the k, the number of fourier terms, that minimizes AIC. In all of my attempts at ARIMA with regressors, it always produces an error when k>3 (or i>3 if we're talking about inside my loop). The error being Error in solve.default(res$hessian * n.used, A) : system is computationally singular: reciprocal condition number = 1.39139e-34. Simply setting k=3 gives some decent results for my test set but for the next year it doesn't appear to adequately catch the steep drops around the end of the year and is much smoother than imagined as evidenced in this forecast:AutoArima with Covariates (k=3)

I assume this general smoothness is due to the small number of fourier pairs. Is there an oversight in my code in that I'm just royally screwing up the procedure provided by Dr. Hyndman? Or is there a theoretical issue that I'm unknowingly running into by trying to find more than 3 pairs of fourier terms for the multiple seasons I'm attempting to account for? Is there a better way to include the multiple seasonalities and dummy variables?

Any help in getting these covariates into the arima model with an appropriate number of fourier terms would be appreciated. If not, I'd at least like to know whether or not what I'm attempting is possible in general with larger number of fourier pairs.

Improve this question
$\endgroup$
6
  • $\begingroup$ What country is the data from ? It appears that you are not taking into account different kinds of patterns around known events or any level/step shifts that might be present or any anomalies that might be distorting patterns thus effecting parameter estimation . Additionally there may be particular days of the month that might be important or long-weekends or particular weeks in the month. Major flaws in using fourier procedures/scripts is the assumption that effects ( like daily effects ) are invariant over time whereas it is quite possible that these effects have changed over time. $\endgroup$
    – IrishStat
    Oct 9, 2015 at 0:00
  • 1
    $\begingroup$ The error message seems to indicate that you have perfect multicollinearity. For example, if you had monthly and quartely Fourier terms, there clearly would be perfect multicollinearity as three months match one quarter. Maybe something analogous is happening in your case, too. $\endgroup$
    – Richard Hardy
    Oct 9, 2015 at 9:31
  • $\begingroup$ You are not looking for outliers. With no outlier identification and adjustment, you won't be able to clearly identify the patterns. $\endgroup$
    – Tom Reilly
    Oct 9, 2015 at 12:03
  • $\begingroup$ @IrishStat: This is US data. Items like long weekends and particular weeks/months are additional dummy variables that I would like to add in once I have the fourier side taken care of. $\endgroup$
    – JRW
    Oct 9, 2015 at 12:44
  • $\begingroup$ @RichardHardy: That was my first instinct even though I wasn't completely familiar with the theory going on behind the scenes. Since I am only using weekly and annual (365.25 days) fourier terms I struggle to see how I could be ending up with perfect multicollinearity. Since this isn't the first time I've ran into this issue I also question how my other data could be perfectly multicollinear as well. $\endgroup$
    – JRW
    Oct 9, 2015 at 12:49

3 Answers 3

Reset to default
8
$\begingroup$

You're hitting the wall because you're exhausting limitations of the first fourier transform fourier(1:n,i,m1). As RandomDude correctly pointed out above, # of transforms i should be less than half period (m1).

However, if, with your code, you run 2 cycles -- one for i, and another for j, where j would be # of transforms for the second seasonality cycle fourier(1:n,j,m3), you would still have a lot of room for model improvement.

This is what I've got from your data, even without dummies, only based on AR, MA, and data seasonality:

library(forecast)
y <- msts(ts, c(7,365)) # multiseasonal ts
fit <- auto.arima(y, seasonal=F, xreg=fourier(y, K=c(3,30)))
fit_f <- forecast(fit, xreg= fourierf(y, K=c(3,30), 180), 180)
plot(fit_f)

enter image description here

I suspect the performance will even improve when holidays are added.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Indeed. After catching my mistake on the looping after @RandomDude's answer, everything is looking good. Even if I had known about the period/2 maximum k value, I should have been keen enough to realize my mistake of subjecting both fourier series to the same k value. Perhaps I should start drinking more coffee... $\endgroup$
    – JRW
    Oct 14, 2015 at 20:25
  • $\begingroup$ Did you see my explanation abt n/2 threshold below? BTW, thanks for the beautiful dataset! What do the data measure? Sales? Load? $\endgroup$
    – Sergey Bushmanov
    Oct 14, 2015 at 20:40
  • $\begingroup$ Yes I did, and I appreciate it. I now regret not doing some more digging around before troubling you guys with the question. In regards to the dataset I'm glad you find it as pretty as I do. Without going into too much detail the measure is dealing with load volume. $\endgroup$
    – JRW
    Oct 14, 2015 at 21:37
  • $\begingroup$ How do you guys make sure to actually find the best k for a fourier term? When you only use a simple for-loop and break-up as soon as the actual AICc is worse than the one last round you may end up with a fourier term that is not optimal. At least i had this issue several times already when i used the method Rob J Hyndman suggested... $\endgroup$
    – RandomDude
    Oct 19, 2015 at 11:42
  • 1
    $\begingroup$ @RandomDude Do you mean search path has local minima, i.e. after an increase in $AIC_c$, if you proceed further, there is a chance to find a model with even a lesser $AIC_c$ ? $\endgroup$
    – Sergey Bushmanov
    Oct 19, 2015 at 12:19
3
$\begingroup$

Is there a reason why you are not using the fourier() function in the forecast package? When you try to build a fourier term of a seasonal time series object your K must be smaller than period/2. Otherwise you get an error:

fourier(ts(test, frequency=7),4) #3 works, 4+ doesn't
Error in ...fourier(x, K, 1:length(x)) : 
  K must be not be greater than period/2

Quote from ?fourier()

fourier(x,K)

When x is a ts object, the value of K should be an integer and specifies the number of sine and cosine terms to return. Thus, the matrix returned has 2*K columns.

I don't have a theoretical explanation + i don't have enough reputation to write a comment under your post (answer was the only option). Hope i could still help you somehow!

Improve this answer
$\endgroup$
4
  • $\begingroup$ To answer the question about why I'm not using the fourier() function in the forecast package, I didn't know it existed. I've noticed that Dr. Hyndman mentions throughout his blog about the function being in the package however it is not in the reference manual. Being the case, I simply thought he was referencing the fourier function he defined at robjhyndman.com/hyndsight/longseasonality. $\endgroup$
    – JRW
    Oct 12, 2015 at 12:47
  • $\begingroup$ The issue with K not being greater than period/2 most certainly is what I'm running into (especially since my loop was trying to calculate K based on both periods simultaneously). Now if only someone smarter than me could explain what the theoretical violation is when period/2. $\endgroup$
    – JRW
    Oct 12, 2015 at 12:51
  • 2
    $\begingroup$ Perhaps I'm not smarter than you, but, from spectral analysis any ts, regardless stochastic or deterministic, with seasonal patterns or without, can be fit with n/2 sin and cos pairs for even n and with (n-1)/2 for odd n, where n is your claimed period. This is why singularity warnings when you try to feed your model more than n/2 predictors. $\endgroup$
    – Sergey Bushmanov
    Oct 13, 2015 at 19:54
  • 2
    $\begingroup$ @JRW this blog post is quite new and deals with the method, maybe you can learn something there link $\endgroup$
    – RandomDude
    Oct 13, 2015 at 22:25
1
$\begingroup$

Optimization of Fourier pairs based on AICc values. This is for yearly and monthly seasonality on data without weekends. The ranges 0:10 and 1:20 should be changed accordingly for different seasonal periods. Or increased for a broader search. 

msts_test <- msts( test , seasonal.periods = c(21.66,260))


my_aic_df <- matrix(ncol = 10 , nrow = 20)


for(i in 1:10){ 

   for(j in 1:20){ 

   fn <- fourier( msts_test , K=  c(i , j) )

   FourierFit <- auto.arima( msts_test , seasonal=FALSE,  xreg=fn )

   my_aic_df[(j),(i+1)] <- FourierFit$aicc

   }
}

 which(my_aic_df == min(my_aic_df), arr.ind = TRUE)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.