5
$\begingroup$

I'm trying to create an Arima model and forecast it ahead the next 20 hours using the code and data below.

When I look at the median of df$tri for each hour and broken down by day of the week, each weekday seems to have a distinct 24 hour pattern.

So I thought I would try adding dummy variables for the day of the week to my model as xreg predictors.

When I plot Acast$mean, I'm just getting a flat line. So I was wondering if there was something incorrect about the way I'm creating and using the dummy variables for day of the week.

Code:

##BoxCox

Tlambda <- BoxCox.lambda(df$Tri)

##Partitioning Time Series

EndTrain<-336
ValStart<-EndTrain+1
ValEnd<-ValStart+20

tsTrain <-df$Tri[1:EndTrain]
tsValidation<-df$Tri[ValStart:ValEnd]


##Weekday variables
Copydf<-df
Copydf$Weekday<-as.factor(weekdays(as.Date(Copydf$DateTime, "%Y-%m-%d")))
Weekdays <- Copydf$Weekday[1:nrow(Copydf)]
xreg1 <- model.matrix(~as.factor(Weekdays)+0)[, 1:7]
colnames(xreg1) <- c("Friday", "Monday", "Saturday", "Sunday","Thursday","Tuesday","Wednesday")

xreg1Train<-xreg1[1:EndTrain,]
xreg1Val<-xreg1[ValStart:ValEnd,]

##Checking effect of Weekday
Arima.fit <- auto.arima(tsTrain, lambda = Tlambda, xreg=xreg1Train, stepwise=FALSE, approximation = FALSE )

##Forecast

Acast<-forecast(Arima.fit,xreg=xreg1Val, h=20)

Acast$mean

Data:

    dput(df$Tri[1:336])
c(11, 14, 17, 5, 5, 5.5, 8, NA, 5.5, 6.5, 8.5, 4, 5, 9, 10, 11, 
7, 6, 7, 7, 5, 6, 9, 9, 6.5, 9, 3.5, 2, 15, 2.5, 17, 5, 5.5, 
7, 6, 3.5, 6, 9.5, 5, 7, 4, 5, 4, 9.5, 3.5, 5, 4, 4, 9, 4.5, 
6, 10, NA, 9.5, 15, 9, 5.5, 7.5, 12, 17.5, 19, 7, 14, 17, 3.5, 
6, 15, 11, 10.5, 11, 13, 9.5, 9, 7, 4, 6, 15, 5, 18, 5, 6, 19, 
19, 6, 7, 7.5, 7.5, 7, 6.5, 9, 10, 5.5, 5, 7.5, 5, 4, 10, 7, 
5, 12, 6, NA, 4, 2, 5, 7.5, 11, 13, 7, 8, 7.5, 5.5, 7.5, 15, 
7, 4.5, 9, 3, 4, 6, 17.5, 11, 7, 6, 7, 4.5, 4, 4, 5, 10, 14, 
7, 7, 4, 7.5, 11, 6, 11, 7.5, 15, 23.5, 8, 12, 5, 9, 10, 4, 9, 
6, 8.5, 7.5, 6, 5, 8, 6, 5.5, 8, 11, 10.5, 4, 6, 7, 10, 11.5, 
11.5, 3, 4, 16, 3, 2, 2, 8, 4.5, 7, 4, 8, 11, 6.5, 7.5, 17, 6, 
6.5, 9, 12, 17, 10, 5, 5, 9, 3, 8.5, 11, 4.5, 7, 16, 11, 14, 
6.5, 15, 8.5, 7, 6.5, 11, 2, 2, 13.5, 4, 2, 16, 11.5, 3.5, 9, 
16.5, 2.5, 4.5, 8.5, 5, 6, 7.5, 9.5, NA, 9.5, 8, 2.5, 4, 12, 
13, 10, 4, 6, 16, 16, 13, 8, 12, 19, 19, 5.5, 8, 6.5, NA, NA, 
NA, 15, 12, NA, 6, 11, 8, 4, 2, 3, 4, 10, 7, 5, 4.5, 4, 5, 11.5, 
12, 10.5, 4.5, 3, 4, 7, 15.5, 9.5, NA, 9.5, 12, 13.5, 10, 10, 
13, 6, 8.5, 15, 16.5, 9.5, 14, 9, 9.5, 11, 15, 14, 5.5, 6, 14, 
16, 9.5, 23, NA, 19, 12, 5, 11, 16, 8, 11, 9, 13, 6, 7, 3, 5.5, 
7.5, 19, 6.5, 5.5, 4.5, 7, 8, 7, 10, 11, 13, NA, 12, 1.5, 7, 
7, 12, 8, 6, 9, 15, 9, 3, 5, 11, 11, 8, 6, 3, 7.5)

    dput(df$DateTime[1:336])
c("2015-01-01 00:00", "2015-01-01 01:00", "2015-01-01 02:00", 
"2015-01-01 03:00", "2015-01-01 04:00", "2015-01-01 05:00", "2015-01-01 06:00", 
"2015-01-01 07:00", "2015-01-01 08:00", "2015-01-01 09:00", "2015-01-01 10:00", 
"2015-01-01 11:00", "2015-01-01 12:00", "2015-01-01 13:00", "2015-01-01 14:00", 
"2015-01-01 15:00", "2015-01-01 16:00", "2015-01-01 17:00", "2015-01-01 18:00", 
"2015-01-01 19:00", "2015-01-01 20:00", "2015-01-01 21:00", "2015-01-01 22:00", 
"2015-01-01 23:00", "2015-01-02 00:00", "2015-01-02 01:00", "2015-01-02 02:00", 
"2015-01-02 03:00", "2015-01-02 04:00", "2015-01-02 05:00", "2015-01-02 06:00", 
"2015-01-02 07:00", "2015-01-02 08:00", "2015-01-02 09:00", "2015-01-02 10:00", 
"2015-01-02 11:00", "2015-01-02 12:00", "2015-01-02 13:00", "2015-01-02 14:00", 
"2015-01-02 15:00", "2015-01-02 16:00", "2015-01-02 17:00", "2015-01-02 18:00", 
"2015-01-02 19:00", "2015-01-02 20:00", "2015-01-02 21:00", "2015-01-02 22:00", 
"2015-01-02 23:00", "2015-01-03 00:00", "2015-01-03 01:00", "2015-01-03 02:00", 
"2015-01-03 03:00", "2015-01-03 04:00", "2015-01-03 05:00", "2015-01-03 06:00", 
"2015-01-03 07:00", "2015-01-03 08:00", "2015-01-03 09:00", "2015-01-03 10:00", 
"2015-01-03 11:00", "2015-01-03 12:00", "2015-01-03 13:00", "2015-01-03 14:00", 
"2015-01-03 15:00", "2015-01-03 16:00", "2015-01-03 17:00", "2015-01-03 18:00", 
"2015-01-03 19:00", "2015-01-03 20:00", "2015-01-03 21:00", "2015-01-03 22:00", 
"2015-01-03 23:00", "2015-01-04 00:00", "2015-01-04 01:00", "2015-01-04 02:00", 
"2015-01-04 03:00", "2015-01-04 04:00", "2015-01-04 05:00", "2015-01-04 06:00", 
"2015-01-04 07:00", "2015-01-04 08:00", "2015-01-04 09:00", "2015-01-04 10:00", 
"2015-01-04 11:00", "2015-01-04 12:00", "2015-01-04 13:00", "2015-01-04 14:00", 
"2015-01-04 15:00", "2015-01-04 16:00", "2015-01-04 17:00", "2015-01-04 18:00", 
"2015-01-04 19:00", "2015-01-04 20:00", "2015-01-04 21:00", "2015-01-04 22:00", 
"2015-01-04 23:00", "2015-01-05 00:00", "2015-01-05 01:00", "2015-01-05 02:00", 
"2015-01-05 03:00", "2015-01-05 04:00", "2015-01-05 05:00", "2015-01-05 06:00", 
"2015-01-05 07:00", "2015-01-05 08:00", "2015-01-05 09:00", "2015-01-05 10:00", 
"2015-01-05 11:00", "2015-01-05 12:00", "2015-01-05 13:00", "2015-01-05 14:00", 
"2015-01-05 15:00", "2015-01-05 16:00", "2015-01-05 17:00", "2015-01-05 18:00", 
"2015-01-05 19:00", "2015-01-05 20:00", "2015-01-05 21:00", "2015-01-05 22:00", 
"2015-01-05 23:00", "2015-01-06 00:00", "2015-01-06 01:00", "2015-01-06 02:00", 
"2015-01-06 03:00", "2015-01-06 04:00", "2015-01-06 05:00", "2015-01-06 06:00", 
"2015-01-06 07:00", "2015-01-06 08:00", "2015-01-06 09:00", "2015-01-06 10:00", 
"2015-01-06 11:00", "2015-01-06 12:00", "2015-01-06 13:00", "2015-01-06 14:00", 
"2015-01-06 15:00", "2015-01-06 16:00", "2015-01-06 17:00", "2015-01-06 18:00", 
"2015-01-06 19:00", "2015-01-06 20:00", "2015-01-06 21:00", "2015-01-06 22:00", 
"2015-01-06 23:00", "2015-01-07 00:00", "2015-01-07 01:00", "2015-01-07 02:00", 
"2015-01-07 03:00", "2015-01-07 04:00", "2015-01-07 05:00", "2015-01-07 06:00", 
"2015-01-07 07:00", "2015-01-07 08:00", "2015-01-07 09:00", "2015-01-07 10:00", 
"2015-01-07 11:00", "2015-01-07 12:00", "2015-01-07 13:00", "2015-01-07 14:00", 
"2015-01-07 15:00", "2015-01-07 16:00", "2015-01-07 17:00", "2015-01-07 18:00", 
"2015-01-07 19:00", "2015-01-07 20:00", "2015-01-07 21:00", "2015-01-07 22:00", 
"2015-01-07 23:00", "2015-01-08 00:00", "2015-01-08 01:00", "2015-01-08 02:00", 
"2015-01-08 03:00", "2015-01-08 04:00", "2015-01-08 05:00", "2015-01-08 06:00", 
"2015-01-08 07:00", "2015-01-08 08:00", "2015-01-08 09:00", "2015-01-08 10:00", 
"2015-01-08 11:00", "2015-01-08 12:00", "2015-01-08 13:00", "2015-01-08 14:00", 
"2015-01-08 15:00", "2015-01-08 16:00", "2015-01-08 17:00", "2015-01-08 18:00", 
"2015-01-08 19:00", "2015-01-08 20:00", "2015-01-08 21:00", "2015-01-08 22:00", 
"2015-01-08 23:00", "2015-01-09 00:00", "2015-01-09 01:00", "2015-01-09 02:00", 
"2015-01-09 03:00", "2015-01-09 04:00", "2015-01-09 05:00", "2015-01-09 06:00", 
"2015-01-09 07:00", "2015-01-09 08:00", "2015-01-09 09:00", "2015-01-09 10:00", 
"2015-01-09 11:00", "2015-01-09 12:00", "2015-01-09 13:00", "2015-01-09 14:00", 
"2015-01-09 15:00", "2015-01-09 16:00", "2015-01-09 17:00", "2015-01-09 18:00", 
"2015-01-09 19:00", "2015-01-09 20:00", "2015-01-09 21:00", "2015-01-09 22:00", 
"2015-01-09 23:00", "2015-01-10 00:00", "2015-01-10 01:00", "2015-01-10 02:00", 
"2015-01-10 03:00", "2015-01-10 04:00", "2015-01-10 05:00", "2015-01-10 06:00", 
"2015-01-10 07:00", "2015-01-10 08:00", "2015-01-10 09:00", "2015-01-10 10:00", 
"2015-01-10 11:00", "2015-01-10 12:00", "2015-01-10 13:00", "2015-01-10 14:00", 
"2015-01-10 15:00", "2015-01-10 16:00", "2015-01-10 17:00", "2015-01-10 18:00", 
"2015-01-10 19:00", "2015-01-10 20:00", "2015-01-10 21:00", "2015-01-10 22:00", 
"2015-01-10 23:00", "2015-01-11 00:00", "2015-01-11 01:00", "2015-01-11 02:00", 
"2015-01-11 03:00", "2015-01-11 04:00", "2015-01-11 05:00", "2015-01-11 06:00", 
"2015-01-11 07:00", "2015-01-11 08:00", "2015-01-11 09:00", "2015-01-11 10:00", 
"2015-01-11 11:00", "2015-01-11 12:00", "2015-01-11 13:00", "2015-01-11 14:00", 
"2015-01-11 15:00", "2015-01-11 16:00", "2015-01-11 17:00", "2015-01-11 18:00", 
"2015-01-11 19:00", "2015-01-11 20:00", "2015-01-11 21:00", "2015-01-11 22:00", 
"2015-01-11 23:00", "2015-01-12 00:00", "2015-01-12 01:00", "2015-01-12 02:00", 
"2015-01-12 03:00", "2015-01-12 04:00", "2015-01-12 05:00", "2015-01-12 06:00", 
"2015-01-12 07:00", "2015-01-12 08:00", "2015-01-12 09:00", "2015-01-12 10:00", 
"2015-01-12 11:00", "2015-01-12 12:00", "2015-01-12 13:00", "2015-01-12 14:00", 
"2015-01-12 15:00", "2015-01-12 16:00", "2015-01-12 17:00", "2015-01-12 18:00", 
"2015-01-12 19:00", "2015-01-12 20:00", "2015-01-12 21:00", "2015-01-12 22:00", 
"2015-01-12 23:00", "2015-01-13 00:00", "2015-01-13 01:00", "2015-01-13 02:00", 
"2015-01-13 03:00", "2015-01-13 04:00", "2015-01-13 05:00", "2015-01-13 06:00", 
"2015-01-13 07:00", "2015-01-13 08:00", "2015-01-13 09:00", "2015-01-13 10:00", 
"2015-01-13 11:00", "2015-01-13 12:00", "2015-01-13 13:00", "2015-01-13 14:00", 
"2015-01-13 15:00", "2015-01-13 16:00", "2015-01-13 17:00", "2015-01-13 18:00", 
"2015-01-13 19:00", "2015-01-13 20:00", "2015-01-13 21:00", "2015-01-13 22:00", 
"2015-01-13 23:00", "2015-01-14 00:00", "2015-01-14 01:00", "2015-01-14 02:00", 
"2015-01-14 03:00", "2015-01-14 04:00", "2015-01-14 05:00", "2015-01-14 06:00", 
"2015-01-14 07:00", "2015-01-14 08:00", "2015-01-14 09:00", "2015-01-14 10:00", 
"2015-01-14 11:00", "2015-01-14 12:00", "2015-01-14 13:00", "2015-01-14 14:00", 
"2015-01-14 15:00", "2015-01-14 16:00", "2015-01-14 17:00", "2015-01-14 18:00", 
"2015-01-14 19:00", "2015-01-14 20:00", "2015-01-14 21:00", "2015-01-14 22:00", 
"2015-01-14 23:00")
$\endgroup$
7
$\begingroup$

To start with we will explore different ways that repeating patterns can appear in time series data and how we can model those patterns. This may be over kill for the question, however I do think that this answer will help you think about what is happening in the models and design better experiments to model your data going forward.

Simple daily seasonality

To begin with, lets think about a daily seasonal ARIMA model. This type of model is looking for some type of pattern that repeats where we see the same thing every day. A time series that this type of model might work well for might look like this:

set.seed(0)
# create a pattern that repeats every 24 hours
daily <- sinpi(seq(1, 3 - 1/24, length.out = 24))
# create a time series out of that repeating pattern 2 weeks long
dailyseas <- ts(rep(daily, 14), frequency = 24)
# add some noise to the data to make it interesting to fit
dailyseas <- dailyseas + rnorm(length(dailyseas), sd = 0.2)
plot(dailyseas, main = "Daily Repeating Pattern")

daily seasonal

We can then fit a seasonal AR model to the data and get a pretty good forecast for the series. Since we know that the pattern is pretty stable over time, I will use two seasonal lags rather than just one because this will allow the model to smooth out any noise in the data better.

> (dsmodel <- arima(dailyseas, seasonal = c(2L, 0L, 0L)))

Call:
arima(x = dailyseas, seasonal = c(2L, 0L, 0L))

Coefficients:
        sar1    sar2  intercept
      0.4254  0.5396     0.0094
s.e.  0.0455  0.0465     0.1368

sigma^2 estimated as 0.05491:  log likelihood = -20.56,  aic = 49.13

daily seasonal forecast

This is just about perfect, the intercept is near zero, which it should be, and the sum of the coefficients on the seasonal lags is close to 1, meaning the forecast is about an average of them (eg. we predict the value tomorrow at noon is about the average of today and yesterday at noon). It is important to note we let the ARIMA model know how often the pattern repeated itself by making a ts object and setting the frequency = 24. Alternatively, we could have used a vector for the series and set seasonal = list(order = c(1L, 0L, 0L), period = 24).

This works well when we have a simple repeating pattern, but what if we have a day of week effect.

Daily Seasonality with a Day of Week effect

A day of week effect is an consistant impact on the underlying series we see based on the day of the week. We can add a day of week effect to our data using:

# Create an adjustment for day of week: we will leave Monday as 0
# so later it is easy to see the change in other days relative to 
# Monday.
(doweffect <- c(0 , sample(c(-3, -2, -1, 1, 2, 3))))
# add our day of week effect to the original series
dowseas <- dailyseas + rep(doweffect, 2, each = 24)
plot(dowseas)

seasonal with day of week effect

We handle this new pattern in our data in one of two ways 1) adding external regressors to our original ARIMA model or 2) thinking of the weekly repeating pattern in the data as the new seasonality of the data. In an ARIMA model with external regressors, we are looking for some sort of ARIMA type pattern that is "thrown off" by some amount by the things quantified by the external regressors. In the model with weekly seasonality, we are looking for the interaction of the daily and weekly pattern and ignoring that the daily pattern is still present in each day of the week. Below we create the original seasonal model as well as the 2 variants.

# Creates a model matrix to indicate the day of week for values in
# our time series.  Note the model matrix does not have a column to 
# indicate Monday.  The purpose of the model matrix is to allow the 
# model to include the impact of a value relative to some baseline, 
# usually the first factor level; in this case, Monday.  We also 
# remove the intercept term since intercept is already part of the 
# ARIMA models.
> dowreg <- model.matrix(~as.factor(rep(1:7, 2, each = 24)))[, -1]
> colnames(dowreg) <- c("Tues", "Weds", "Thurs", "Fri", "Sat", "Sun")
# Creates a first model where the day of week is ignored
> (dowsmodel <- arima(dowseas, seasonal = c(2L, 0L, 0L)))

Call:
arima(x = dowseas, seasonal = c(2L, 0L, 0L))

Coefficients:
        sar1     sar2  intercept
      0.1085  -0.1550     0.0170
s.e.  0.0588   0.0612     0.1123

sigma^2 estimated as 4.455:  log likelihood = -728.46,  aic = 1464.93

Hmm, that doesn't look great, our intercept is near zero and our seasonal coefficients are nearly canceling each other out. Lets look at our external regressor model:

# Creates a model where the day of week effect is accounted for using 
# an external regressor
> (dowxreg <- arima(dowseas, seasonal = c(2L, 0L, 0L), xreg = dowreg))

Call:
arima(x = dowseas, seasonal = c(2L, 0L, 0L), xreg = dowreg)

Coefficients:
        sar1    sar2  intercept    Tues    Weds    Thurs      Fri      Sat     Sun
      0.4301  0.5351    -0.0080  1.0376  2.0036  -0.9759  -1.9841  -3.0051  3.0363
s.e.  0.0459  0.0468     0.1389  0.0394  0.0360   0.0424   0.0427   0.0413  0.0438

sigma^2 estimated as 0.05457:  log likelihood = -19.51,  aic = 59.01

This is much better, the two seasonal lags (sar1 and sar2) are basically taking an average again like they did in our other model, and the external regressors are adjusting the days by the correct amount (1 for Tues, 2 for Weds, -1 for Thurs, -2 for Fri, -3 for Sat and 3 for Sun). How about our weekly seasonal model:

# Creates a model where the day of week effect is accounted for by 
# increasing the seasonality to weekly rather than daily.  This time 
# we can only use 1 seasonal lag because our data only don't have 
# enough seasonal periods at the weekly frequency.
> (dowlongs <- arima(dowseas, seasonal = list(order = c(1L, 0L, 0L), 
+                                             period = 24*7)))

Call:
arima(x = dowseas, seasonal = list(order = c(1L, 0L, 0L), period = 24 * 7))

Coefficients:
      sar1  intercept
         1      0.005
s.e.   NaN   1400.371

sigma^2 estimated as 0.03082:  log likelihood = -8.64,  aic = 23.27

Once again, this looks good, it is predicting next Monday at noon should be the same as last Monday at noon, which makes sense for this data. Lets look at how this turns out in the forecasts:

par(mfrow = c(1,1))
plot(forecast(dowsmodel, 48), PI = FALSE, xlim = c(8, 17), 
     main = "Forecasts from Various Seasonal AR models")
lines(forecast(dowxreg, 48, xreg = dowreg[1:48, ])$mean, col = "red", lwd = 2)
lines(forecast(dowlongs, 48)$mean, col = "green", lwd = 2)
abline(v = 8, lty = 2)
legend("topleft", bty = "n",
       c("Seasonal AR", "Seasonal AR with DoW", "Long Seasonal AR"),
       fill = c("blue", "red", "green"))

forecasts DoW

While the original mode that worked so well before falls apart, we see that both of the other approaches work well for this new data. The model with external regressors is able to find the daily pattern that is occurring once it takes into account the effect of the day of week. The weekly seasonal model is missing the daily pattern, but is able to overcome that by seeing the larger pattern that repeats every week.

Distinct pattern for each day of the week

Now we are finally going to get into what you claim to see in your data; a daily pattern which is different for each day of the week. We can make a series with this property as follows:

# Creates a time series where each day of week has a unique pattern
dailysig <- c(daily, 
              sort(daily),
              sort(daily, decreasing = TRUE),
              abs(daily),
              exp(daily),
              log(daily + 2),
              cos(daily))
# Creates a two week long version of this series with a noise component
diffdaily <- ts(rep(dailysig, 2), frequency = 24)
diffdaily <- diffdaily + rnorm(length(diffdaily), sd = 0.2)
plot(diffdaily, main = "Unique Pattern for Each Day of Week")

weekly seasonal

We see that Monday's pattern doesn't really look like Tuesday's or Wednesday's etc. Lets examine what happens if we try and make each of our 3 types of models for this data set.

# Makes the same three models as last time
> (dowsmodel <- arima(diffdaily, seasonal = c(2L, 0L, 0L)))

Call:
arima(x = diffdaily, seasonal = c(2L, 0L, 0L))

Coefficients:
        sar1     sar2  intercept
      0.2488  -0.2669     0.4651
s.e.  0.0533   0.0537     0.0390

sigma^2 estimated as 0.509:  log likelihood = -365.57,  aic = 739.14

> (dowxreg <- arima(diffdaily, seasonal = c(2L, 0L, 0L), xreg = dowreg))

Call:
arima(x = diffdaily, seasonal = c(2L, 0L, 0L), xreg = dowreg)

Coefficients:
        sar1     sar2  intercept     Tues     Weds   Thurs     Fri     Sat     Sun
      0.0756  -0.3189     0.0224  -0.0670  -0.0105  0.5791  1.1814  0.6315  0.7423
s.e.  0.0531   0.0532     0.0875   0.1205   0.1419  0.1241  0.1199  0.1325  0.1234

sigma^2 estimated as 0.3413:  log likelihood = -298.76,  aic = 617.52

> (dowlongs <- arima(diffdaily, seasonal = list(order = c(1L, 0L, 0L), period = 24*7)))

Call:
arima(x = diffdaily, seasonal = list(order = c(1L, 0L, 0L), period = 24 * 7))

Coefficients:
        sar1  intercept
      0.9254     0.4591
s.e.  0.0111     0.0575

sigma^2 estimated as 0.08293:  log likelihood = -221.5,  aic = 448.99

> plot(forecast(dowsmodel, 48), PI = FALSE, xlim = c(8, 17), 
+      main = "Forecasts from Various Seasonal AR models for Different DoW Effects")
> lines(forecast(dowxreg, 48, xreg = dowreg[1:48, ])$mean, col = "red", lwd = 2)
> lines(forecast(dowlongs, 48)$mean, col = "green", lwd = 2)
> abline(v = 8, lty = 2)
> legend("topleft", bty = "n",
+        c("Seasonal AR", "Seasonal AR with DoW", "Long Seasonal AR"),
+        fill = c("blue", "red", "green"))

weekly seasonal forecast

Results fell apart for the first two models this time. Since there isn't an underlying daily pattern any more, the model with external regressors was only able to see that certain days of the week are higher or lower on average, but the pattern from hour to hour was missed. The weekly seasonal model however was still about to see the weekly repeat and make a reasonable model.

Your Data

Now that we have seen the importance of seasonality in our models, lets see what happens if we try running auto.arima again, but this time making your data a seasonal time series.

> tsTrain <- ts(tsTrain, frequency = 24)
> (dowsmodel <- auto.arima(tsTrain))
Series: tsTrain 
ARIMA(0,0,0)(1,0,0)[24] with non-zero mean 

Coefficients:
        sar1    mean
      0.0508  8.4899
s.e.  0.0579  0.2452

sigma^2 estimated as 17.31:  log likelihood=-928.96
AIC=1863.91   AICc=1863.98   BIC=1875.36
> (dowxreg <- auto.arima(tsTrain, xreg = dowreg))
Series: tsTrain 
Regression with ARIMA(0,0,0)(1,0,0)[24] errors 

Coefficients:
        sar1  intercept     Tues    Weds   Thurs     Fri     Sat     Sun
      0.0841     7.7401  -0.4165  2.3504  0.0211  1.1671  1.8975  0.3029
s.e.  0.0587     0.5991   0.8074  0.8496  0.8617  0.8520  0.8461  0.8288

sigma^2 estimated as 16.63:  log likelihood=-919.55
AIC=1857.1   AICc=1857.65   BIC=1891.45
> (dowlongs <- auto.arima(ts(tsTrain, frequency = 24*7)))
Series: ts(tsTrain, frequency = 24 * 7) 
ARIMA(0,0,2) with non-zero mean 

Coefficients:
         ma1      ma2    mean
      0.2433  -0.0171  8.5118
s.e.  0.0583   0.0506  0.2782

sigma^2 estimated as 16.51:  log likelihood=-921.06
AIC=1850.12   AICc=1850.24   BIC=1865.39
> plot(forecast(dowsmodel, 24), PI = FALSE, xlim = c(8, 16), 
+      main = "Forecasts from Various Seasonal AR models for Different DoW Effects")
> lines(forecast(dowxreg, 24, xreg = dowreg[1:24, ])$mean, col = "red", lwd = 2)
> lines(ts(forecast(dowlongs, 24)$mean, start = 15, frequency = 24), 
+       col = "green", lwd = 2)
> abline(v = 8, lty = 2)
> legend("topleft", bty = "n",
+        c("Seasonal AR", "Seasonal AR with DoW", "Long Seasonal AR"),
+        fill = c("blue", "red", "green"))

OP forecast

The "each weekday seems to have a distinct 24 hour pattern" doesn't seem to be happening as seen by the trouble fitting a the weekly seasonal model, but there does seem to be a daily seasonality the models are picking up on since. Personally, I would trust the plain seasonal model (no external regressors) the most since it is less prone to over fitting than the one with external regressors, but that is your call. In general, you might feel disappointed since the forecasts don't look much like your data. This is because there is a lot of noise in your data that the model still can't account for.

Conclusions

  1. A seasonal model will allow your data to find repeating patterns in your data.
  2. Adding external regressors to your model can allow the model to find the underlying pattern when the pattern is obscured by another influence.
  3. If every day of the week has a different pattern, that is a weekly seasonality, not a day of week effect.
$\endgroup$
  • $\begingroup$ Why not add both daily effects and hourly effects and allow for unusual values to be detected and remedied with pulse indicators being added to the model where needed ? $\endgroup$ – IrishStat Aug 31 '17 at 17:15
  • $\begingroup$ @IrishStat Do you mean an ARIMA model with two different seasonalities? If so, that does sound appealing, but I am not aware of any implementation for such a model in the languages I am familiar with. You could of course use a different model like tbats which is based on exponential smoothing, but the OP was requesting an ARIMA based model. $\endgroup$ – Barker Aug 31 '17 at 20:37
  • $\begingroup$ i meant a deterministic model with 1 + 6 + 23 predictors and then identifying/incorpoarating an approriate arima model and potential outliers/level shifts/seasonal pulses/time trends AND arma. IN this way the OP gets his arma strtcture bur it is integrated into needed deterministic structure, there is no double seasonality in ARIMA. $\endgroup$ – IrishStat Aug 31 '17 at 21:53
  • $\begingroup$ @IrishStat I would hesitate to put a model with so many predictors on this size dataset. In particular, if the pattern really does show a differently pattern every day of the week (which I am not convince it does), you would need 24*7 - 1 predictors to account for the fact that "2pm" has a different pattern on Monday than on Tuesday. If you want to account for multiple seasonalities, I would be more inclined to try including Fourier terms to avoid parameter space explosion. $\endgroup$ – Barker Aug 31 '17 at 22:58
  • $\begingroup$ If you use 6 day-of-the-week and 23 hour of the day that's 29 = 1 for the constant. You can search/try to see what interaction series are significant )which essentially forms a minimally sufficient set . Then when you add any needed ARIMA structure you are home along with intervention detection to ferret out unusual values. I do this routinely with great results. . $\endgroup$ – IrishStat Sep 1 '17 at 0:35
0
$\begingroup$

Hyndman's docs say the xreg vector needs to have the same number of rows as the time series. In your code, where defining 'Weekdays' you are missing a comma before the closing square bracket.

If this external regressor approach doesn't work I'd try fitting a seasonal arima model with m=7 manually.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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