Добавление сезонных манекенов при составлении прогноза разложения STL и модели ARIMA

Question

В настоящее время я прогнозирую структуру трафика, анализируя данные из API Google Distance Matrix. Престижность этого поста в блоге ( https://petolau.github.io/Forecast-electricity-consumption-with-similar-day-approach-in-R/ Similar-day-approach-in-R/) и Роба Дж. Хиндмана за пакет прогнозов, который я использую.

Я определил набор траекторий и путей, которые представляют собой набор траекторий. Каждая из этих траекторий будет прогнозироваться для каждого дня недели с понедельника по воскресенье отдельно, путем объединения разложения STL для прогнозирования сезонной модели и модели ARIMA для прогнозирования остатка. Я думаю, что результаты довольно хороши, но я не могу проверить это, так как данные для Германии и в настоящее время есть школьные и праздничные дни.

Я использую stl и forecast функции. Поскольку этот подход был самым многообещающим из-за всего, что я тестировал, я бы очень хотел придерживаться его. В настоящее время мои прогнозы в час пик намного выше, чем прогнозы Google. Есть ли способ добавить сезонные манекены в качестве регрессоров к этой модели, чтобы она не переоценивала трафик? Большое спасибо.

Набор данных: https://drive.google.com/file/d/0BzptouUUdw_jbHZENW1WVGt4Y2s/view?usp=sharing

Код:

library(data.table)
library(lubridate)
library(ggplot2)
library(forecast)
library(openxlsx)
library(imputeTS)

#load excel file
wb <- loadWorkbook("Data.xlsx")

#write sheets into a list of data frames
name <- names(wb)
P1 <- c("t1","t2","t3","t8","t9","t10")
P2 <- c("t1","t2","t4","t5","t6","t7","t9","t10")
P3 <- c("t1","t2","t4","t11","t12","t13")
P <- list("P1" = P1, "P2"=P2, "P3"=P3)

#all values
df <- lapply(name, function(x) readWorkbook(wb, sheet = x, startRow = 1, colNames = TRUE, 
                  rowNames = FALSE, detectDates = TRUE, skipEmptyRows = TRUE, skipEmptyCols = TRUE, 
                  rows = NULL, cols = NULL, check.names = FALSE, namedRegion = NULL, 
                  na.strings = "NA", fillMergedCells = FALSE))

#assign names to columns
names(df) <- name

#merge DFs into one list
DT <- rbindlist(df)

#convert columns
DT[, datetime := as.POSIXct( DT$datetime*24*3600 + as.POSIXct("1899-12-30 00:00", 
                                                tz="Etc/GMT-1"), tz="Etc/GMT-1")]
DT[, date := as.Date(DT$datetime, tz="Etc/GMT-1", format= "%d-%m-%Y")]
DT[, value := gsub("[^0-9]", "", DT$value)]
DT[, value := as.numeric(DT$value)]
DT[, value := na.locf(DT$value, option = "locf", na.remaining = "rev")]

#check how many segments contain full data
count_ID <- DT[, .N, date]
full <- count_ID[N == max(N), .(date)]
nrow(full) # number of extracted IDs
period <- 144

#add weekdays
DT[, week := weekdays(datetime)]
unique(DT[, week])

#extract relevant vectors
n_weekdays <- unique(DT[, week])
n_date <- unique(DT[, date])

unique(DT[, ID]) # all road segments



## create an overview for the sections
date <- data.table(DT[ID %in% "t1", datetime])
sP <- list(date)
for (i in 1:length(P))
{
  sP[[i+1]] <- lapply(P[[i]], function(x) data.table(DT[ID %in% x, value]))
}

sP1 <- data.table(as.data.table(sP[[1]]), as.data.table(sP[[2]]))
sP1[ ,sum := rowSums(.SD), .SDcols = 2:7]
sP1[, date := as.Date(sP1$V1, tz="Etc/GMT-1", format= "%d-%m-%Y")]
sP2 <- data.table(as.data.table(sP[[1]]), as.data.table(sP[[3]]))
sP2[ ,sum := rowSums(.SD), .SDcols = 2:9]
sP2[, date := as.Date(sP2$V1, tz="Etc/GMT-1", format= "%d-%m-%Y")]
sP3 <- data.table(as.data.table(sP[[1]]), as.data.table(sP[[4]]))
sP3[ ,sum := rowSums(.SD), .SDcols = 2:7]
sP3[, date := as.Date(sP3$V1, tz="Etc/GMT-1", format= "%d-%m-%Y")]

## STL + ARIMA
stlARIMAPred <- function(Y, period = 144){ #dummyperiod
  ts_Y <- ts(Y, start = 0, freq = period)
  dekom <- stl(ts_Y, s.window = "periodic", robust = TRUE)
  arima <- forecast(dekom, h = period, method = "arima")
  return(as.vector(arima$mean))
}


#create function that predicts forecast for a week
predictWeek <- function(data, set_of_date, FUN, train_win = 0){

  for_mon <- FUN(data[(week == n_weekdays[4] & date %in% set_of_date), value])
  for_tue <- FUN(data[(week == n_weekdays[5] & date %in% set_of_date), value])
  for_wed <- FUN(data[(week == n_weekdays[6] & date %in% set_of_date), value])
  for_thu <- FUN(data[(week == n_weekdays[7] & date %in% set_of_date), value])
  for_fri <- FUN(data[(week == n_weekdays[1] & date %in% set_of_date), value])
  for_sat <- FUN(data[(week == n_weekdays[2] & date %in% set_of_date), value])
  for_sun <- FUN(data[(week == n_weekdays[3] & date %in% set_of_date), value])

  return(c(for_mon, for_tue, for_wed, for_thu, for_fri, for_sat, for_sun))
}


## create forecast for road segments and write them into data table

Arimas <- lapply(name, function (x) predictWeek(DT[ID %in% x], n_date[1:24], stlARIMAPred) ) 
Arimas <- as.data.table(Arimas)
names(Arimas) <- name

# create list with different paths

ArimaP <- list()
for (i in 1:length(P)){
  ArimaP[[i]] <- lapply(P[[i]], function(x) data.table(Arimas[,get(x)]))
}


sAP1 <- as.data.table(ArimaP[[1]])
sAP1[ ,sum := rowSums(.SD), .SDcols = 1:6]
sAP2 <- as.data.table(ArimaP[[2]])
sAP2[ ,sum := rowSums(.SD), .SDcols = 1:8]
sAP3 <- as.data.table(ArimaP[[3]])
sAP3[ ,sum := rowSums(.SD), .SDcols = 1:6]

Информация о сессии:

R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)