Smoothing methods
Сорокина Валерия
05 декабря 2018
1. Загрузка данных
Ряд данных (TRP_M_PASS_DIRI) содержит месячные показатели индекса пассажирооборота транспорта общего пользования
Пассажирооборот транспорта общего пользования - объем работы транспорта по перевозкам пассажиров. Единицей измерения является пассажиро-километр. Исчисляется суммированием произведений количества пассажиров каждой перевозки на расстояние перевозки в километрах. Перевозки грузов, грузооборот и пассажирооборот транспорта определяются на основании сводных итогов, представляемых МПС России, Минтрансом России, ОАО “Газпром”, АК “Транс-нефтепродукт” и данных Госкомстата России по централизованной государственной статистиче-ской отчетности по видам транспорта и оценке объемов транспортной работы, выполненной пред-приятиями, относящимися к субъектам малого предпринимательства и предпринимателями (физическими лицами), занимающимися грузовыми и пассажирскими перевозками.
tr <- sophisthse("TRP_M_PASS_DIRI")
sophisthse_metadata(tr)## # A tibble: 1 x 7
## tsname unit fullname methodology source comment freq
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 TRP_M_… 1997.… Индекс пасса… "Пассажирооборо… Федерал… Пересчет и… 122. Визуализация ряда
dygraph(tr, main = "Пассажирооборот транспорта общего пользования",
xlab = "Год",
ylab = "пассажиро-километр") %>%
dyLegend(show = "follow") %>%
dyRangeSelector()Наблюдается ярко выраженная сезонность и прослеживается тренд. Повышенный спрос на пассажирские перевозки наблюдается в летние месяцы, в то время как пониженный спрос - в зимние. Обусловлено это может быть тем, что в летние месяцы население более склонно в отпускам, что приводит к повышению спроса на транспортные услуги. Более того, летом многие выезжают на дачи или пикники, для чего также зачастую используется общественный транспорт.
tr %>%
decompose(type = 'multiplicative') %>%
autoplot() +
labs(title = 'Декомпозиция мультипликативного ряда')
3. Влияние кризиса
Первый спад наблюдается в 1998 году, на который приходится тяжелейший кризис в стране.
Второй спад произошел в 2009 году, когда наблюдался мировой финансовый упадок. Активная фаза кризиса завершилась в декабре 2009, далее спрос восстанавливается и наблюдается небольшой рост
Третий спад приходится на 2015 - 2016 гг., что может объясняться кризисом, развивавшимся в 2014 и 2015 годах в РФ.
Изменений в сезонности конкретно в кризисные года не наблюдалось. Однако, на протяжении всего рассматриваемого периода можно заметить колебания сезонных коэффициентов при построении по методу Винтерса:
tr_ob <- window(tr, end = c(2017,10))
tr_ob_winters_m <- hw(tr_ob, h = 12, seasonal = 'multiplicative')
autoplot(cbind('Сезонность' = tr_ob_winters_m$model$states[, 's1']), facets = T) +
labs(title = "Пассажирооборот транспорта общего пользования: сезонность", y = NULL, x = NULL, color = 'Обозначения')
4. Выбор модели прогнозирования
Так как в ряде данных наблюдается выраженные тренд и сезонность, наиболее подходящей моделью для прогнозирования станет Метод Винтерса.
Метод Винтерса позволяет учитывать как скорость изменения ряда, так и сезонность. Существует две модификации метода Винтерса, которые отличаются подходом к учету сезонности: мультипликативная и аддитивная модели.
5. Органичение глубины исторического периода
На первый взгляд амплитуда колебаний непостоянна, а увеличивается с течением времени. Можно предполодить, что необходима мультипликативная модель. Для более корректного анализа следует оставить весь доступный исторический период.
6. Ex-Post прогнозирование
Для выбора лучшей модели прогнозирования на предварительном этапе сравниваются аддитивная и мультипликативная модели метода Винтерса и модели с разными значениями констант.
В качестве тестового периода используются последние 12 месяцев от имеющихся данных.
tr_ob <- window(tr, end = c(2017,10))
tr_test <- window(tr, start = c(2017,11))
dygraph(cbind(tr_ob, tr_test), main = "Пассажирооборот транспорта общего пользования") %>%
dySeries("tr_ob", label = "Обучающие") %>%
dySeries("tr_test", label = "Тестовые") %>%
dyRangeSelector() %>%
dyLegend(show = "follow")Мультипликативная модель Винтерса
tr_ob_winters_m <- hw(tr_ob, h = 12, seasonal = 'multiplicative') # Константы сглаживания и начальное значение уровня и тренда подбираются автоматически
# График прогноза
autoplot(tr_ob_winters_m) +
autolayer(tr_test, series = "Тестовый ряд") +
labs(title = "Пассажирооборот транспорта общего пользования",
fill = 'Доверительный\nинтервал',
x = NULL, y = NULL)
Коэффициенты модели:
tr_ob_winters_m$model$par %>% round(3)## alpha beta gamma l b s0 s1 s2 s3
## 0.230 0.000 0.528 119.916 0.200 0.976 0.905 0.962 0.992
## s4 s5 s6 s7 s8 s9 s10
## 1.216 1.278 1.184 1.012 0.910 0.911 0.798Аддитивная модель Винтерса
tr_ob_winters_a <- hw(tr_ob, h = 12, seasonal = 'additive') # Константы сглаживания и начальное значение уровня и тренда подбираются автоматически
# График прогноза
autoplot(tr_ob_winters_a) +
autolayer(tr_test, series = "Тестовый ряд") +
labs(title = "Пассажирооборот транспорта общего пользования",
fill = 'Доверительный\nинтервал',
x = NULL, y = NULL)
Коэффициенты модели:
tr_ob_winters_a$model$par %>% round(3)## alpha beta gamma l b s0 s1 s2 s3
## 0.381 0.000 0.595 117.306 0.146 -11.872 -9.583 -1.902 4.631
## s4 s5 s6 s7 s8 s9 s10
## 32.995 39.336 24.040 0.681 -11.096 -19.745 -26.205Изменение константы
Можно заметить, что мультипликативная модель Винтерса визуально более близка к факту, нежели аддитивная, но все же не полностью совпадает. Коэффициент beta = 0, то есть модель никак не учитывает тренд. Поэтому введем коэффициент вручную. Коэффициент при этом возьмем приближенный к нулю, так как резкого скачка в данных нет.
tr_ob_winters_beta <- hw(tr_ob, alpha = 0.23, beta = 0.05, h = 12, seasonal = 'multiplicative')
# График прогноза
autoplot(tr_ob_winters_beta) +
autolayer(tr_test, series = "Тестовый ряд") +
labs(title = "Пассажирооборот транспорта общего пользования",fill = 'Доверительный\nинтервал',
x = NULL, y = NULL)
Сравнение моделей
Визуальное сравнение
Исходя из визуального сравнения, наиболее приближенным вариантом является мультипликативная модель с поправкой beta.
window(tr, start = 2017) %>%
autoplot(tr_ob, series = 'Обучающий') +
autolayer(tr_test, series = 'Тестовый') +
autolayer(tr_ob_winters_m, series = 'Мультипликативная', PI = FALSE) +
autolayer(tr_ob_winters_a, series = 'Аддитивная', PI = FALSE) +
autolayer(tr_ob_winters_beta, series = 'Поправка beta', PI = FALSE) +
labs(title = "Пассажирооборот транспорта общего пользования", color = NULL, x = NULL, y = NULL)
Сравнение ошибок
rbind(accuracy(tr_ob_winters_m),
accuracy(tr_ob_winters_a),
accuracy(tr_ob_winters_beta)) %>%
as_tibble() %>%
round(4) %>%
mutate(`Метод` = c('Мультипликативная модель',
'Аддитивная модель',
'Поправка beta')) %>%
select(`Метод`, ME, MAPE, MPE) %>%
arrange(ME)## # A tibble: 3 x 4
## Метод ME MAPE MPE
## <chr> <dbl> <dbl> <dbl>
## 1 Мультипликативная модель -0.231 2.72 -0.356
## 2 Аддитивная модель -0.0218 2.74 -0.174
## 3 Поправка beta 0.0819 2.97 -0.0684Отличия в показателях ошибки незначительные. Наилучшее согласие с данными демонстрирует мультипликативная модель с поправкой beta.
7. Прогноз на 2 года
tr_winters_m <- hw(tr, beta = 0.05, h = 24, seasonal = 'multiplicative')
# График прогноза
autoplot(tr_winters_m) +
labs(title = "Пассажирооборот транспорта общего пользования",
fill = 'Доверительный\nинтервал',
x = NULL, y = NULL)
8. Оценка стандартной ошибки прогноза
При помощи метода скользящего контроля оценим стандартную ошибку прогноза на горизонте 1, 3, 6, 12 и 24 месяцев:
RMSE <- function(errors) {
errors^2 %>% mean(na.rm = T) %>% sqrt() }
tr %>% tsCV(hw, h = 1) %>% RMSE() %>% round(2)## [1] 6.26tr %>% tsCV(hw, h = 3) %>% RMSE() %>% round(2)## [1] 9.95tr %>% tsCV(hw, h = 6) %>% RMSE() %>% round(2)## [1] 12.95tr %>% tsCV(hw, h = 12) %>% RMSE() %>% round(2)## [1] 16.98tr %>% tsCV(hw, h = 24) %>% RMSE() %>% round(2)## [1] 26.97Можно заметить: чем больше горизонт прогнозиования, тем больше ошибка.
Источники