Чтобы оценить, насколько эффективно работают методы прогнозирования на определенных данных, применяются показатели ошибки прогноза. В этом блокноте рассмотрены несколько групп таких показателей - абсолютные, относительные, сравнительные, которые с различных сторон характеризуют ошибку прогноза, т.е. отклонение прогнозируемых и фактических значений.
Наиболее часто для оценки эффективности прогнозирования используется относительный (процентный) показатель ошибки - MAPE, однако у данного показателя есть ряд недостатков, которые проявляются в случае, когда временной ряд содержит нулевые или очень малые наблюдения. Даже когда нет проблемы нулевых значений, качество прогнозирования сложно оценить лишь с использованием MAPE, поскольку для каждого конкретного объекта прогнозирования действует свой набор факторов (новизна продукта, маркетинговая активность, количество и тип клиентов, степень агрегирования), который определяет пределы точности прогнозирования. Поэтому нельзя ориентироваться на усредненные, “отраслевые” значения ошибки прогнозирования и необходимо сравнивать подходы к прогнозированию в одинаковых условиях. С этой целью будут рассмотрены показатели, основанные на сопоставлении ошибки прогноза применяемого для прогнозирования метода с показателями ошибки простых, “наивных” методов прогнозирования на тех же данных.
Также будут рассмотрены методы оценки качества моделей прогнозирования, основанные на разделении данных на “обучающий” и “тестовый” периоды, а также скользящем контроле.
Материал составлен на основе раздела 3.1 и 3.4 книги: Hyndman R, Athanasopoulos G. Forecasting Principles and Practice.
library(tidyverse) # визуализация и трансформация данных
library(forecast) # анализ временных рядов и прогнозирование
library(fpp) # Примеры временных рядов
library(sophisthse) # Загрузка временных рядов из базы Sophist
Вначале рассмотрим несколько простых методов прогнозирования, которые могут использоваться как основа для сравнения при бенчмаркинге различных методов прогнозирования.
В этом методе в качестве прогноза на будущие периоды берется среднее значение ряда, вычисленное на всем историческом периоде:
\[ \hat{y}_{T+h} = \frac{1}{T} \sum_{t=1}^{T} y_t \]
Метод также может применяться и для неупорядоченных данных (берется среднее по выборке).
В R для реализации этого метода используется функция: forecast::meanf(y, h), где y - временной ряд, а h - горизонт прогнозирования.
Наивный прогноз - это простое повторение последнего наблюдения в качестве прогноза на будущие периоды. В R метод реализован двумя функциями: forecast::naive(y, h) и forecast::rwf(y, h).
Этот метод использует тот же принцип, однако повторяется не последнее наблюдение, а последнее наблюдение для одноименного периода. Например, прогноз на декабрь 2018 года, сделанный в 2017 году, будет равен значению за декабрь 2017 года.
В R этот метод реализован функцией forecast::snaive(y, h)
Дрейф - это медленное изменение уровня ряда. В методе дрейфа вычисляется средний прирост уровня ряда по всем смежным периодам. Затем это значение используется для прогнозирования будущих значений.
\[ \hat{y}_{T+h} = y_T + \frac{h}{T-1} \sum_{t=2}^{T} (y_t - y_{t-1})\]
Можно показать, что данное выражение эквивалентно:
\[ \hat{y}_{T+h} = y_T + h \frac{y_T - y_1}{T-1} \]
Т.е. в методе дрейфа прогнозом является продолжение прямой линии, проведенной через первое и последнее наблюдения.
В R метод дрейфа реализован функцией: forecast::rwf(y, h, drift=TRUE).
# Отбор данных
beer2 <- ausbeer %>%
window(start = 1992, end = c(2005, 4))
beer_label <- 'Поквартальный прогноз производства пива в Австралии'
beer_unit <- 'млн. л'
# Прогнозирование
beer_m <- meanf(beer2, h = 11) # по среднему значению ряда
beer_n <- naive(beer2, h = 11) # наивный прогноз
beer_sn <- snaive(beer2, h = 11) # сезонный наивный прогноз
beer_d <- rwf(beer2, h = 11, drift = T) # метод дрейфа
Для визуализации прогнозов можно воспользоваться функциями autoplot() и autolayer():
# Сравнение нескольких прогнозов
autoplot(beer2, series = 'История') +
autolayer(beer_m, PI = FALSE, series = 'Среднее ряда') +
autolayer(beer_n, PI = FALSE, series = 'Наивный прогноз') +
autolayer(beer_sn, PI = FALSE, series = 'С. наивный прогноз') +
autolayer(beer_d, PI = FALSE, series = 'Метод дрейфа') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
Пакет forecast рассчитывает и доверительные интервалы для прогноза, которые легко получить с в виде графика или таблицы. В предыдущем примере мы отключили вывод доверительных интервалов для прогноза с помощью параметра PI = FALSE, чтобы не загромождать график.
# График
autoplot(beer_sn) +
labs(title = beer_label, y = beer_unit,
x = NULL, color = NULL)
# Таблица
beer_sn
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2006 Q1 416 393.9027 438.0973 382.2051 449.7949
## 2006 Q2 403 380.9027 425.0973 369.2051 436.7949
## 2006 Q3 408 385.9027 430.0973 374.2051 441.7949
## 2006 Q4 482 459.9027 504.0973 448.2051 515.7949
## 2007 Q1 416 384.7497 447.2503 368.2068 463.7932
## 2007 Q2 403 371.7497 434.2503 355.2068 450.7932
## 2007 Q3 408 376.7497 439.2503 360.2068 455.7932
## 2007 Q4 482 450.7497 513.2503 434.2068 529.7932
## 2008 Q1 416 377.7264 454.2736 357.4655 474.5345
## 2008 Q2 403 364.7264 441.2736 344.4655 461.5345
## 2008 Q3 408 369.7264 446.2736 349.4655 466.5345
Доверительные интервалы позволяют оценить степень надежности прогноза. Чем шире интервал, тем выше неопределенность.
Объект, который возвращают функции для прогнозирования, имеет следующую структуру:
names(beer_d)
## [1] "method" "model" "lambda" "x" "fitted"
## [6] "residuals" "series" "mean" "level" "lower"
## [11] "upper"
Наиболее полезными компонентами являются mean (точечный прогноз), fitted (прогноз в историческом периоде) и residuals - остатки прогноза. С назначением остальных компонентов можно познакомиться в справке: ?forecast.
Для выделения прогноза в историческом периоде и остатков удобно использовать функции fitted() и residuals() соответственно.
# Подогнанные значения
autoplot(beer2, series = 'История') +
autolayer(fitted(beer_sn), series = 'Прогноз') +
labs(title = beer_label, y = beer_unit,
x = NULL, color = NULL)
## Warning: Removed 4 rows containing missing values (geom_path).
# Остатки
autoplot(residuals(beer_sn)) +
labs(title = 'Ошибки прогноза объема производства пива',
y = beer_unit, x = NULL, color = NULL)
Рассчитаем прогноз показателя из ЕАЭСД (sophist): Ввод в действие жилых домов, млн.кв.метров (CONSTR_Y_NAT)
constr <- sophisthse('CONSTR_Y_NAT')
#sophisthse_metadata(constr)
constr_label <- 'Ежегодный ввод в действие жилых домов (CONSTR_Y_NAT)'
constr_unit <- 'млн.кв.м'
# Оставляем данные с 2005 года
constr <- constr %>%
window(start = 2005)
# Визуализация
autoplot(constr) +
labs(title = constr_label, y = constr_unit, x = NULL)
Прогнозирование ввода новых жилых домов с помощью рассмотренных методов.
constr_m <- meanf(constr, h = 10)
constr_n <- naive(constr, h = 10)
constr_d <- rwf(constr, h = 10, drift = T)
autoplot(constr, series = 'История') +
autolayer(constr_m, PI = FALSE, series = 'Среднее ряда') +
autolayer(constr_n, PI = FALSE, series = 'Наивный прогноз') +
autolayer(constr_d, PI = FALSE, series = 'Метод дрейфа') +
labs(title = constr_label, x = NULL, color = NULL, y = constr_unit)
Ошибка прогноза в периоде \(t\) - это разность фактических и прогнозных значений:
\[ e_t = y_t -\hat{y}_t \]
Ошибку прогноза можно вычислить для каждого прошедшего периода. Для получения общей картины о точности прогнозирования, применяются показатели ошибки прогнозирования, которые вычисляются с помощью агрегирования ошибок за все периоды.
Все показатели ошибки прогнозирования можно разделить на две группы:
абсолютные - зависят от масштаба пронозируемой величины (scale-dependent) и, как правило, измеряются в тех же единицах, что и эта величина;
относительные - представляют собой безразмерные отношения и не зависят от масштаба (scale-independent). Показатели этой группы удобно использовать для сравнения точности прогнозирования для разных объектов.
Еще одним критерием для классификации показателей может служить их назначение. Одни показатели характеризуют разброс прогнозных и фактических значений, другие - систематическую ошибку, т.е. смещение прогноза.
Наиболее популярными абсолютными показателями являются средняя абсолютная ошибка (Mean Absolute Error, MAE) и стандартная ошибка (Root Mean Squared Error, RMSE):
\[ MAE = mean(|e_t|) \]
\[ RMSE = \sqrt{mean(e_t^2)} \]
Оба показателя позволяют сравнивать точность прогнозирования при работе с одним и тем же рядом. MAE легче вычисляется и интерпретируется, а RMSE применяется для построение доверительных интервалов и расчета страхового запаса.
Кроме того, к абсолютным относится еще один показатель - средняя ошибка (Mean Error, ME), который позволяет оценить систематическую ошибку (смещение) прогноза:
\[ ME = mean(e_t) \]
Наиболее популярным относительным показателем является средняя абсолютная ошибка в процентах (Mean Absolute Percentage Error, MAPE). Этот показатель вычисляется путем усреднения относительных ошибок прогноза:
\[ p_t = 100 \frac{e_t} {y_t} \]
\[ MAPE = mean(|p_t|) \]
Для оценки систематической ошибки прогноза применяется показатель средняя ошибка в процентах (Mean Percentage Error, MPE):
\[ MPE = mean(p_t) \]
Относительные ошибки удобны тем, что с их помощью можно сравнивать точность прогнозирования для разных временных рядов, отличающихся масштабом величин.
Недостатком этих показателей является невозможность их вычисления, если встречаются нулевые значения ряда (например, при прогнозировании спроса на малооборачиваемые или сезонные товары). Величина MAPE получается очень большой, когда фактические значения спроса \(y_t < 1\).
MAPE и MPE нельзя применять, когда у показателя нет абсолютного нуля (например, при прогнозировании температуры в шкале Цельсия), поскольку для таких данных не имеет смысла само понятие отношения (деления).
Другим способом нормирования ошибок с целью получения относительных показателей является сопоставление их с абсолютной ошибкой наивного метода (scaling). Для несезонных рядов используется деление на среднюю абсолютную ошибку наивного прогноза, для сезонных - используется сезонный наивный прогноз.
Нормированная ошибка (scaled error) для несезонного ряда:
\[ q_j = \frac{e_j} {\frac{1}{T - 1} \sum_{t=2}^T |y_t - y_{t-1}|} \]
Нормированная ошибка для сезонного ряда:
\[ q_j = \frac{e_j} {\frac{1}{T - m} \sum_{t=m+1}^T |y_t - y_{t-m}|} \]
Для неупорядочынных (перекрестных) данных нормированная ошибка вычисляется относительно стандартного отклонения прогнозируемой величины:
\[ q_j = \frac{e_j} {\frac{1}{N} \sum_{i=1}^N |y_i - \bar{y}|} \]
Средняя абсолютная нормированная ошибка (Mean Absolute Scaled Error) вычисляется путем усреднения нормированных ошибок для всех периодов:
\[ MASE = mean(|q_j|) \]
Для вычисления показателей ошибки в R используется функция forecast::accuracy().
Ряд сезонный, поэтому наилучший результат среди простых методов должен давать сезонный наивный прогноз.
# Метод усреднения
autoplot(beer2, series = 'История') +
autolayer(fitted(beer_m), series = 'Среднее ряда') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
# Наивный прогноз
autoplot(beer2, series = 'История') +
autolayer(fitted(beer_n), series = 'Наивный прогноз') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
## Warning: Removed 1 rows containing missing values (geom_path).
# Сезонный наивный прогноз
autoplot(beer2, series = 'История') +
autolayer(fitted(beer_sn), series = 'С. наивный прогноз') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
## Warning: Removed 4 rows containing missing values (geom_path).
# Метод дрейфа
autoplot(beer2, series = 'История') +
autolayer(fitted(beer_d), series = 'Метод дрейфа') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
## Warning: Removed 1 rows containing missing values (geom_path).
Показатели ошибки прогноза автоматически рассчитываются функцией forecast::accuracy():
accuracy(beer_m)
## ME RMSE MAE MPE MAPE MASE
## Training set 8.121418e-15 44.1763 35.91135 -0.9510944 7.995509 2.444228
## ACF1
## Training set -0.1256697
Для удобства сравнения прогнозов, объединим результаты по всем прогнозам в одну таблицу:
rbind(accuracy(beer_m),
accuracy(beer_n),
accuracy(beer_sn),
accuracy(beer_d)) %>%
as_tibble() %>%
round(2) %>%
mutate(`Метод` = c('Среднее', 'Наивный', 'С. наивный', 'Дрейф')) %>%
select(`Метод`, MASE, MAPE, everything()) %>%
arrange(MASE)
## # A tibble: 4 x 8
## Метод MASE MAPE ME RMSE MAE MPE ACF1
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 С. наивный 1 3.4 -1.85 17.2 14.7 -0.48 -0.34
## 2 Среднее 2.44 8 0 44.2 35.9 -0.95 -0.13
## 3 Наивный 3.77 12.3 0.71 66.6 55.4 -0.9 -0.25
## 4 Дрейф 3.78 12.3 0 66.6 55.5 -1.06 -0.25
Показатель ошибки MASE для сезонного наивного прогноза равен 1, поскольку для сезонного ряда (когда задана частота) именно этот прогноз используется как основа для сравнения. Усреднение в ~2.5 раза хуже по этому показателю, чем сезонный наивный прогноз.
Примечание: по умолчанию выводится слишком большое количество значащих цифр. Для упрощения интерпретации показателей можно округлить все значения в таблице с помощью функции round(). Обратите внимание, что это надо делать до того, как в таблицу будет добавлен текстовый столбец с названиями методов.
При выборе метода прогнозирования аналитика интересует не то, как эта модель объясняла прошлое, а то, как она предсказывает будущее. Поэтому показатели ошибки, вычисленные на историческом периоде могут вводить в заблуждение. Можно построить такую модель, которая “предсказывает” уже известное прошлое идеально точно.
Поэтому распространенным приемом для выбора модели является разделение ряда на два периода - обучающий (training set) и тестовый (test set). При этом модель строится на обучающем периоде, а показатели ошибки вычисляются на тестовом. Этот способ называется прогнозирование Ex-Post.
Прогнозирование Ex-Post
Обычно в качестве тестового множества выбирается порядка 20% наблюдений, однако это ориентировочная величина, которую можно скорректировать исходя из объема имеющихся данных и горизонта прогнозирования. Рекомендуется делать тестовый период как минимум такой же длины, как и горизонт прогнозирования.
В примере о прогнозировании пива модели строились на подмножестве ряда ausbeer Для проверки доступно еще 11 ежеквартальных наблюдений.
autoplot(ausbeer, series = 'Полный ряд') +
autolayer(beer2, series = 'Сокращенный ряд')
Прогнозы в тестовом периоде можно сравнить с фактическими данными визуально.
# Сравнение нескольких прогнозов
autoplot(window(ausbeer, start = 2004), series = 'Факт') +
autolayer(beer_m, PI = FALSE, series = 'Среднее ряда') +
autolayer(beer_n, PI = FALSE, series = 'Наивный прогноз') +
autolayer(beer_sn, PI = FALSE, series = 'С. наивный прогноз') +
autolayer(beer_d, PI = FALSE, series = 'Метод дрейфа') +
labs(title = beer_label, x = NULL, color = NULL, y = beer_unit)
Функция accuracy() автоматически вычисляет показатели ошибки на тестовом периоде, если предоставить данные тестового периода, или просто весь ряд целиком.
# Метод усреднения
accuracy(beer_m, ausbeer)
## ME RMSE MAE MPE MAPE MASE
## Training set 8.121418e-15 44.17630 35.91135 -0.9510944 7.995509 2.444228
## Test set -1.718344e+01 38.01454 33.77760 -4.7345524 8.169955 2.298999
## ACF1 Theil's U
## Training set -0.12566970 NA
## Test set -0.08286364 0.7901651
# Наивный прогноз
accuracy(beer_n, ausbeer)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.7090909 66.60207 55.43636 -0.8987351 12.26632 3.773156
## Test set -62.2727273 70.90647 63.90909 -15.5431822 15.87645 4.349833
## ACF1 Theil's U
## Training set -0.25475212 NA
## Test set -0.08286364 1.428524
# Сезонный наивный прогноз
accuracy(beer_sn, ausbeer)
## ME RMSE MAE MPE MAPE MASE
## Training set -1.846154 17.24261 14.69231 -0.4803931 3.401224 1.0000000
## Test set -2.545455 12.96849 11.27273 -0.7530978 2.729847 0.7672537
## ACF1 Theil's U
## Training set -0.3408329 NA
## Test set -0.1786912 0.22573
# Метод дрейфа
accuracy(beer_d, ausbeer)
## ME RMSE MAE MPE MAPE MASE
## Training set 1.240308e-14 66.59830 55.50083 -1.062644 12.29043 3.777543
## Test set -6.652727e+01 74.83196 67.64793 -16.567964 16.79620 4.604310
## ACF1 Theil's U
## Training set -0.25475212 NA
## Test set -0.07101826 1.509823
При использовании тестового периода функция вычисляет еще один показатель - U-статистику Тейла (Theil’s U).
Этот показатель похож на MASE, и в нем используется для нормализации показатель ошибки наивного прогноза. Однако, в отличие от MASE, используются не абсолютные а процентные ошибки. Поэтому данный показатель страдает от тех же недостатков при малых или нулевых значениях ряда, что и показатели MAPE и MPE.
Суммы квадратов процентных ошибок (Sum of Squared Percent Errors, SSPE) для модели и наивного прогноза вычисляются по формулам:
\[ SSPE(модель) = \sum_{t=2}^N \left(\frac{e_t} {y_{t-1}} \right)^2 \]
\[ SSPE(\text{наивный прогноз}) = \sum_{t=1}^N \left( \frac{y_t - y_{t-1}} {y_{t-1}} \right)^2 \]
\[ U= \sqrt{ \frac{SSPE(модель)}{SSPE(\text{наивный прогноз}) }} \]
Интерпретация U-статистики Тейла аналогична интерпретации показателя MASE:
Чтобы сравнить модели на основе ошибки в тестовом периоде, можно собрать все показатели ошибки в одну таблицу, с которой можно работать, используя функции пакета dplyr:
rbind(accuracy(beer_m, ausbeer)['Test set',],
accuracy(beer_n, ausbeer)['Test set',],
accuracy(beer_sn, ausbeer)['Test set',],
accuracy(beer_d, ausbeer)['Test set',]) %>%
as_tibble() %>%
round(2) %>%
mutate(`Метод` = c('Среднее', 'Наивный', 'С. наивный', 'Дрейф')) %>%
select(`Метод`, MASE, MAPE, everything()) %>%
arrange(MASE)
## # A tibble: 4 x 9
## Метод MASE MAPE ME RMSE MAE MPE ACF1 `Theil's U`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 С. наивный 0.77 2.73 -2.55 13.0 11.3 -0.75 -0.18 0.23
## 2 Среднее 2.3 8.17 -17.2 38.0 33.8 -4.73 -0.08 0.79
## 3 Наивный 4.35 15.9 -62.3 70.9 63.9 -15.5 -0.08 1.43
## 4 Дрейф 4.6 16.8 -66.5 74.8 67.6 -16.6 -0.07 1.51
Более сложный и полезный метод оценки качества прогноза основан на скользящем контроле (cross-validation, или rolling forecasting origin). При этом в качестве обучающего периода используются \(n\) начальных наблюдений, а в качестве тестового - одно наблюдение, отстоящее от обучающего периода на величину горизонта прогнозирования (\(n + h\)). Значение \(n\) постепенно увеличивается, насколько позволяют данные. Ошибки прогноза на \(h\) шагов вперед усредняются. Принцип работы метода показан на рисунке:
Скользящий контроль для оценки ошибки прогноза
В это примере горизонт прогноза \(h = 4\). Синие точки - обучающий период, красные - тестовый.
В R скользящий контроль реализован в функции forecast::tsCV(). Эта функция возвращает ошибки прогноза на горизонте \(h\) для каждого периода.
beer3 <- ausbeer %>% window(start = 2003)
tsCV(beer3, meanf, h = 4)
## h=1 h=2 h=3 h=4
## 2003 Q1 -55.00000 -14.000000 55.000000 0.000000
## 2003 Q2 13.50000 82.500000 27.500000 -17.500000
## 2003 Q3 78.00000 23.000000 -22.000000 0.000000
## 2003 Q4 3.50000 -41.500000 -19.500000 22.500000
## 2004 Q1 -42.20000 -20.200000 21.800000 -16.200000
## 2004 Q2 -13.16667 28.833333 -9.166667 -22.166667
## 2004 Q3 30.71429 -7.285714 -20.285714 -15.285714
## 2004 Q4 -11.12500 -24.125000 -19.125000 54.875000
## 2005 Q1 -22.88889 -17.888889 56.111111 12.111111
## 2005 Q2 -15.60000 58.400000 14.400000 -37.600000
## 2005 Q3 59.81818 15.818182 -36.181818 -17.181818
## 2005 Q4 10.83333 -41.166667 -22.166667 63.833333
## 2006 Q1 -42.00000 -23.000000 63.000000 -1.000000
## 2006 Q2 -20.00000 66.000000 2.000000 -42.000000
## 2006 Q3 67.33333 3.333333 -40.666667 -29.666667
## 2006 Q4 -0.87500 -44.875000 -33.875000 45.125000
## 2007 Q1 -44.82353 -33.823529 45.176471 -7.823529
## 2007 Q2 -31.33333 47.666667 -5.333333 -35.333333
## 2007 Q3 49.31579 -3.684211 -33.684211 -13.684211
## 2007 Q4 -6.15000 -36.150000 -16.150000 NA
## 2008 Q1 -35.85714 -15.857143 NA NA
## 2008 Q2 -14.22727 NA NA NA
## 2008 Q3 NA NA NA NA
Для агрегирования ошибок можно использовать любой общепринятый показатель, например стандартную ошибку, RMSE. Для удобства, создадим функцию, вычисляющую этот показатель.
RMSE <- function(errors) {
errors^2 %>% mean(na.rm = T) %>% sqrt()
}
# Усреднение
beer3 %>% tsCV(meanf, h = 4) %>% RMSE() %>% round(1)
## [1] 34.5
# Наивный прогноз
beer3 %>% tsCV(naive, h = 4) %>% RMSE() %>% round(1)
## [1] 48.7
# Сезонный наивный прогноз
beer3 %>% tsCV(snaive, h = 4) %>% RMSE() %>% round(1)
## [1] 23.1
# Метод дрейфа
beer3 %>% tsCV(rwf, drift = TRUE, h = 4) %>% RMSE() %>% round(1)
## [1] 71.2
Исследуем, как зависит стандартная ошибка прогноза от горизонта прогнозирования.
beer_err <- tibble(horizon = 1:12) %>%
mutate(error = map_dbl(horizon,
function(x) tsCV(beer3, naive, h = x) %>% RMSE()))
ggplot(beer_err, aes(horizon, error)) +
geom_line() +
labs(title = 'RMSE прогноза производства пива по наивной модели',
x = 'h', y = beer_unit) +
scale_x_continuous(breaks = beer_err$horizon)