Solved – How to setup xreg argument in auto.arima() in R?

arimatime series

I am working on a small project with one time series which measures the customer visit data (daily). My covariates are a continuous variable Day to measure how many days have been elapsed since the first day of data collection, and some dummy variables, such as whether that day is Christmas, and which day of the week it is, etc.

Part of my data looks like:

Date    Customer_Visit  Weekday Christmas       Day
11/28/11        2535       2        0            1   
11/29/11        3292       3        0            2   
11/30/11        4103       4        0            3   
12/1/11         4541       5        0            4   
12/2/11         6342       6        0            5  
12/3/11         7205       7        0            6   
12/4/11         3872       1        0            7   
12/5/11         3270       2        0            8   
12/6/11         3681       3        0            9

My plan is to use ARIMAX model to fit the data. This can be done in R, with the function auto.arima(). I understand that I have to put my covariates into the xreg argument, but my code for this part always returns an error.

Here is my code:

xreg     <- c(as.factor(modelfitsample$Christmas), as.factor(modelfitsample$Weekday), 
              modelfitsample$Day)
modArima <- auto.arima(ts(modelfitsample$Customer_Visit, freq=7), allowdrift=FALSE, 
                       xreg=xreg)

The error message returned by R is:

Error in model.frame.default(formula = x ~ xreg, drop.unused.levels = TRUE) 
 :variable lengths differ (found for 'xreg')

I learned a lot from How to fit an ARIMAX-model with R? But I am still not very clear how to set up the covariates or dummies in the xreg argument in auto.arima() function.

Best Answer

The main problem is that your xreg is not a matrix. I think the following code does what you want. I've used some artificial data to check that it works.

library(forecast)
# create some artifical data
modelfitsample <- data.frame(Customer_Visit=rpois(49,3000),Weekday=rep(1:7,7),
                             Christmas=c(rep(0,40),1,rep(0,8)),Day=1:49)

# Create matrix of numeric predictors
xreg <- cbind(Weekday=model.matrix(~as.factor(modelfitsample$Weekday)), 
                  Day=modelfitsample$Day,
              Christmas=modelfitsample$Christmas)

# Remove intercept
xreg <- xreg[,-1]

# Rename columns
colnames(xreg) <- c("Mon","Tue","Wed","Thu","Fri","Sat","Day","Christmas")

# Variable to be modelled
visits <- ts(modelfitsample$Customer_Visit, frequency=7)

# Find ARIMAX model
modArima <- auto.arima(visits, xreg=xreg)

Related Solutions

ARIMAX Model – How to Fit an ARIMAX Model Using R for Time Series Analysis

You're going to have a little bit of trouble modeling a series with 2 levels of seasonality using an ARIMA model. Getting this right is going highly dependent on setting things up correctly. Have you considered a simple linear model yet? They're a lot faster and easier to fit than ARIMA models, and if you use dummy variables for your different seasonality levels they are often quite accurate.

I'm assuming you have hourly data, so make sure your TS object is setup with a frequency of 24.
You can model other levels of seasonality using dummy variables. For example, you might want a set of 0/1 dummies representing the month of the year.
Include the dummy variables in the xreg argument, along with any covariates (like temperature).
Fit the model with the arima function in base R. This function can handle ARMAX models through the use of the xreg argument.
Try the Arima and auto.arima functions in the forecast package. auto.arima is nice because it will automatically find good parameters for your arima model. However, it will take FOREVER to fit on your dataset.
Try the tslm function in the arima package, using dummy variables for each level of seasonality. This will fit a lot faster than the Arima model, and may even work better in your situation.
If 4/5/6 don't work, THEN start worrying about transfer functions. You have to crawl before you can walk.
If you are planning to forecast into the future, you will first need to forecast your xreg variables. This is easy for seasonal dummies, but you'll have to think about how to make a good weather forecasts. Maybe use the median of historical data?

Here is an example of how I would approach this:

#Setup a fake time series
set.seed(1)
library(lubridate)
index <- ISOdatetime(2010,1,1,0,0,0)+1:8759*60*60
month <- month(index)
hour <- hour(index)
usage <- 1000+10*rnorm(length(index))-25*(month-6)^2-(hour-12)^2
usage <- ts(usage,frequency=24)

#Create monthly dummies.  Add other xvars to this matrix
xreg <- model.matrix(~as.factor(month))[,2:12]
colnames(xreg) <- c('Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec')

#Fit a model
library(forecast)
model <- Arima(usage, order=c(0,0,0), seasonal=list(order=c(1,0,0), period=24), xreg=xreg)
plot(usage)
lines(fitted(model),col=2)

#Benchmark against other models
model2 <- tslm(usage~as.factor(month)+as.factor(hour))
model3 <- tslm(usage~as.factor(month))
model4 <- rep(mean(usage),length(usage))

#Compare the 4 models
library(plyr) #for rbind.fill
ACC <- rbind.fill(  data.frame(t(accuracy(model))),
                    data.frame(t(accuracy(model2))),
                    data.frame(t(accuracy(model3))),
                    data.frame(t(accuracy(model4,usage)))
                )
ACC <- round(ACC,2)
ACC <- cbind(Type=c('Arima','LM1','Monthly Mean','Mean'),ACC)
ACC[order(ACC$MAE),]

Solved – R Time series forecasting: Having issues selecting fourier pairs for ARIMA with regressors

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)

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

Best Answer

Related Solutions

ARIMAX Model – How to Fit an ARIMAX Model Using R for Time Series Analysis

Solved – R Time series forecasting: Having issues selecting fourier pairs for ARIMA with regressors

Related Question