Использование ARIMA с экзогенными регрессорами для обнаружения выбросов в R

Question

Я хотел бы обнаружить выбросы в данных в реальном времени, которые агрегируются за час.Для этого примера я выбрал почасовые данные о пешеходах из Мельбурна, Австралия ( Объем пешеходов (обновляется ежемесячно) , Система подсчета пешеходов )

IПонимаю, что существует большое количество существующих алгоритмов обнаружения , которые со временем я выучу и использую.В краткосрочной перспективе я хотел бы использовать самый простой подход.Один из таких методов описан @Aksakal в следующей статье об обмене стеками:

Какой алгоритм следует использовать для обнаружения аномалий во временных рядах?

Я думаю,ключ - «неожиданный» квалификатор в вашем графике.Чтобы обнаружить неожиданный , вам нужно иметь представление о том, что ожидается .

. Я бы начал с простой модели временных рядов, такой как AR (p) илиАРМА (р, д).Подгоните его под данные, добавьте сезонность по мере необходимости.Например, ваша модель SAR (1) (24) может быть: $ y_ {t} = c + \ phi y_ {t-1} + \ Phi_ {24} y_ {t-24} + \ Phi_ {25} y_ {t-25} + \ varepsilon_t $, где $ t $ - время в часах.Итак, вы будете предсказывать график на следующий час.Если ошибка прогнозирования $ e_t = y_t- \ hat y_t $ слишком велика, вы выдаете предупреждение.

Когда вы оцениваете модель, вы получите дисперсию $ \ sigma_ \ varepsilon $ ошибки $.\ varepsilon_t $.В зависимости от ваших предположений о распределении, например нормальных, вы можете установить пороговое значение на основе вероятности, например, $ | e_t | <3 \ sigma_ \ varepsilon $ для 99,7% или односторонний $ e_t> 3 \ sigma_ \ varepsilon $.

Количество посетителей, вероятно, довольно постоянное, но супер сезонное.Лучше попробовать сезонные манекены вместо мультипликативной сезонности, тогда вы попробуете ARMAX, где X обозначает экзогенные переменные, которые могут быть чем-то вроде праздничных манекенов, часовых манекенов, манекенов выходных и т. Д.

К сожалению, в посте нет подробностей, поэтому у меня есть несколько вопросов:

Q.1) Как рассчитать / извлечь дисперсию $ \ sigma_ \ varepsilon $ термина ошибки ARIMA $ \ epsilon $из подобранной модели, произведенной auto.arima (данные, xreg = xreg)?Ниже приведен полный пример R, который использует множественную сезонность для определения ежедневной, еженедельной и годовой сезонности.Это не оптимизировано и представлено только в качестве примера реализации, чтобы помочь ответить на вопрос 2.

Я хочу предсказать пороговые значения на весь год (или, по крайней мере, на 30-дневный период)это означает, что h = 24 часа * 30 = 720. По сути, я хочу прогнозировать не среднее значение почасового подсчета пешеходов, а максимальное ожидаемое число пешеходов в час (например, 3σ_ε) для h >> 1 (например,h = 720 часов (30 дней) или даже h = 24 * 365 = 8760 часов (1 год)).

Q.2) Как я могу добиться этого, используя метод, описанный выше?

Пример кода, помогающий решить вышеуказанные вопросы.

library(rwalkr)
library(forecast)
library(tidyverse)
library(tsibble)
library(xts)
library(dygraphs)

pedestrian <- as_tibble(rwalkr::run_melb( year = c(2015:2018) )) 

pedestrian_statelibrary <- pedestrian %>% 
  filter(Sensor == "State Library") %>% 
  left_join(tsibble::holiday_aus(2015:2018, state='VIC'), by=c( 'Date' = 'date' )) %>%
  mutate(holiday = replace_na(holiday, ''),
         Count = ifelse(Count == 0, NA, Count))

# Replace all counts of zero with NA so Box-Cox transform lambda = 0 and constrain output to +ve.
pedestrian_statelibrary_train <- pedestrian_statelibrary %>% filter(Date >= as.Date('2015-05-13'), Date < as.Date('2017-01-01') )
pedestrian_statelibrary_test <- pedestrian_statelibrary %>% filter(Date >= as.Date('2017-01-01') )

# tsbox functions to convert tsibble to tz indirectly. Must be a better way of doing this...
pedestrian_statelibrary_train_zoo <- tsbox::ts_zoo( pedestrian_statelibrary_train %>% select(Date_Time, Count) )
pedestrian_statelibrary_train_ts <-    tsbox::ts_ts(pedestrian_statelibrary_train_zoo)

pedestrian_statelibrary_test_zoo <- tsbox::ts_zoo(    pedestrian_statelibrary_test %>% select(Date_Time, Count) )
pedestrian_statelibrary_test_ts <- tsbox::ts_ts(pedestrian_statelibrary_test_zoo)


## Create external regressors.
xreg_holidays_train <- model.matrix(~as.factor(pedestrian_statelibrary_train$holiday))
xreg_holidays_train <- xreg_holidays_train[,-1]  # remove intercept.
# Remove 1st level from levels()
colnames(xreg_holidays_train) <- levels(as.factor(pedestrian_statelibrary_train$holiday))[-1]

xreg_holidays_test <- model.matrix(~as.factor(pedestrian_statelibrary_test$holiday))
xreg_holidays_test <- xreg_holidays_test[,-1]  # remove intercept.
colnames(xreg_holidays_test) <- levels(as.factor(pedestrian_statelibrary_test$holiday))[-1]

# periods (intervals(samples) per period) for hourly data.
period_day <- 24
period_week <- 24*7
period_year <- 24*365.25

seasonal_periods = c(period_day, period_week, period_year)

pedestrian_statelibrary_train_msts <- msts(pedestrian_statelibrary_train_ts,
                                     start = start(pedestrian_statelibrary_train_ts), 
                                             seasonal.periods = seasonal_periods) 

pedestrian_statelibrary_test_msts <- msts(pedestrian_statelibrary_test_ts, 
                                      start = start(pedestrian_statelibrary_test_ts), 
                                       seasonal.periods = seasonal_periods) 

# set number of Fourier terms per season. Not optimal.
Ks = c(12, 10, 2)

xreg_train <- cbind( seasonality = fourier(pedestrian_statelibrary_train_msts, K = Ks), 
                 holidays = xreg_holidays_train ) 

######################################
## Fit model of exogenous factors and ARIMA as error
######################################
fit <- pedestrian_statelibrary_train_msts %>% 
  auto.arima( xreg = xreg_train,
              seasonal=FALSE,
              stepwise = FALSE,
              parallel = TRUE,
              num.cores = NULL,
              lambda = 0
              ) 

######################################
## Forecast
######################################

fc <- forecast( fit, 
            xreg=cbind( seasonality = fourier(pedestrian_statelibrary_test_msts, K = Ks), 
                        holidays = xreg_holidays_test) 
) 

######################################
## Check residuals and accuracy.
######################################

checkresiduals(fit)

checkresiduals(fc)

accuracy(fc, pedestrian_statelibrary_test_msts)


######################################
## Display fitted model and forecast using interactive dygraph.
######################################

# Plotting `forecast` prediction using `dygraphs`
# https://stackoverflow.com/questions/43624634/plotting-forecast-prediction-using-dygraphs#43668603
as.forecast.ts <- function(forecast_obj){

  training <- forecast_obj$x
  lower <- forecast_obj$lower[,2]
  upper <- forecast_obj$upper[,2]
  point_forecast <- forecast_obj$mean

  cbind(training, lower, upper, point_forecast)
}

fc_ts <- as.forecast.ts(fc)

# Add the time stamps back to ts object.
idx_train <- pedestrian_statelibrary_train %>% ungroup() %>%    select(Date_Time) %>% as.data.frame()
idx_test <- pedestrian_statelibrary_test %>% ungroup() %>% select(Date_Time) %>% as.data.frame()
idx_all <- rbind(idx_train, idx_test)

# Append testing values to fc_ts object, by left joining two xts objects.
test_xts <- as.xts(x = pedestrian_statelibrary_test %>% 
                     dplyr::ungroup() %>%
                     as.data.frame() %>% 
                     dplyr::select( Count ) %>%
                     dplyr::rename( 'testing' = 'Count'), 
                   pedestrian_statelibrary_test$Date_Time)


fc_xts <- as.xts(x = fc_ts %>% 
                   as.data.frame(),
                 idx_all$Date_Time )

fc_xts <- fc_xts %>% xts::merge.xts(test_xts, join='left')

dygraph(data = fc_xts, main = "Pedestrian traffic Forecasting for State Library.") %>% 
  dyRangeSelector %>%
  dySeries(name = "training", label = "Train") %>%
  dySeries(name = 'testing', label = "Test") %>%
  dySeries(name = "point_forecast", label = "Predicted") %>%
  dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
  dyOptions(axisLineColor = "navy", gridLineColor = "grey")