В этом параграфе мы обсудим методы оценки параметров, применимые только к прогнозным моделям. Да, с помощью них можно оценить параметры и регрессий, но ничего кроме более низкой скорости вычисления это не даст.
Итак.
Ошибка прогноза h шагов вперёд
Как мы уже говорили, прогнозирование — это всегда взгляд в будущее. А раз так, решили прогнозисты, то и при оценке моделей стоит использовать критерий, который учитывал бы некоторую ошибку в «будущем».
Здесь стоит сделать небольшое пояснение. Обычно при построении прогнозных моделей (таких, например, как модель экспоненциального сглаживания) на каждом наблюдении делается прогноз на один шаг вперёд \(\hat{y}_{t+1|t}\) и рассчитывается соответствующая одношаговая ошибка \(e_{t+1|t}\). Затем для оценки коэффициентов моделей выбирается какой-нибудь из рассмотренных нами ранее критериев, в него подставляется ряд полученных по всей выборке ошибок, после чего минимизируется. Этот метод используется уже более 60 лет не только в научных исследованиях, но и на практике, при построении, например, моделей экспоненциального сглаживания или авторегрессий. Однако ничто не мешает вместо одношагового делать прогноз на \(h\) шагов вперёд (например, 12), брать эту полученную ошибку и использовать её в каком-нибудь критерии из рассмотренных нами ранее. Например, MSE в этом случае может выглядеть так:
\begin{equation} \label{eq:estim_MSEhsteps}
MSE = \frac{1}{T} \sum_{t=1}^T (y_{t+h} — \hat{y}_{t+h|t})^2 = \frac{1}{T} \sum_{t=1}^T e_{t+h|t}^2 ,
\end{equation}
где \(\hat{y}_{t+h|t}\) — это прогноз из наблюдения \(t\) на \(h\) шагов вперёд, а \(e_{t+h|t}\) — соответствующая прогнозная ошибка.
Использование такого метода позволяет получать более точные прогнозы на выбранный срок (см., например, Cox. Prediction by exponentially weighted moving averages andrelated methods , Clements and Hendry. Multi-step estimation for forecasting), что уже является достаточным аргументом для использования такого метода оценки. Но, конечно же, у него есть и свои минусы.
Во-первых, он более сложен в вычислениях, так как из каждой точки нужно сделать прогноз не на 1, а на \(h\) шагов вперёд.
Во-вторых, если нам нужны прогнозы на несколько наблюдений (например, от 9 до 12), то нам потребуется, следуя вышеизложенной логике, построить несколько моделей — по каждой на срок прогнозирования (4 модели в нашем примере). Это влечёт за собой увеличение времени расчётов. Если стоит задача дать такой прогноз по одному ряду данных, то это увеличение будет несущественным. Если же ставится задача спрогнозировать продажи товара в 10000 разных магазинах, то даже незначительное изменение во времени (например, на 1 секунду) может привести к серьёзным проблемам (увеличение времени расчёта почти на 3 часа)…
В-третьих, из-за прогноза на \(h\) шагов вперёд мы лишаемся \(h\) наблюдений. Например, в выборке из 50 наблюдений при \(h=10\) последнее наблюдение, из которого можно будет дать такой прогноз — это наблюдение номер \(50 — 10 = 40\).
Ну, и в-четвёртых, такой метод оценки является слабо обоснованным со статистической точки зрения. Используя этот метод, мы, фактически, признаём, что либо имеем дело с заведомо неверной, либо неправильно оценённой моделью. К тому же, расчёт различных статистик (например, дисперсии) и вопрос выбора модели в этих условиях оказываются затруднёнными.
Стоит заметить, что использование этого метода при оценке регрессионных моделей (множественная регрессия с экзогенными переменными) обычно не имеет смысла. А всё дело в том, что в регрессиях прогнозные значения зависят не от момента времени \(t\), а от значений переменных, включённых в модель.
Ошибка прогноза от 1 до h шагов вперёд
Похожим по идеи с предыдущим является следующий метод оценки. Вместо того, чтобы записывать прогнозные ошибки только на \(h\) шагов вперёд, можно записать все ошибки от 1 до \(h\), после чего применить к ним любой из понравившихся нам критериев. Например, MSE:
\begin{equation} \label{eq:estim_MSE1tohsteps}
MSE = \frac{1}{h} \sum_{j=1}^h \frac{1}{T} \sum_{t=1}^T (y_{t+j} — \hat{y}_{t+j|t})^2 = \frac{1}{Th} \sum_{j=1}^h \sum_{t=1}^T e_{t+j|t}^2 ,
\end{equation}
В этом случае будет получаться своеобразная средняя во времени оценка, так как в формуле \eqref{eq:estim_MSE1tohsteps} складываются и краткосрочные, и долгосрочные ошибки. В этом как раз и заключается плюс такого метода — оценив модель с помощью \eqref{eq:estim_MSE1tohsteps}, мы получаем прогнозы как на краткосрочную, так и средне- и долгосрочную перспективы. Однако здесь кроется один небольшой изъян — обычно дисперсия долгосрочных ошибок больше дисперсии краткосрочных, так как с ростом горизонта прогнозирования растёт и неопределённость. Это будет приводить к тому, что сумма квадратов ошибок в \eqref{eq:estim_MSE1tohsteps} будет выше для \(j\) близких к \(h\), чем к 1. В связи с этим минимизация такой целевой функции будет приводить в первую очередь к уменьшению долгосрочных ошибок, в ущерб краткосрочным. Однако по скорости вычисления этот критерий оказывается эффективней критерия \eqref{eq:estim_MSEhsteps}, применённого к каждому шагу прогнозирования, причём разница в скорости будет тем выше, чем выше горизонт прогнозирования.
Ряд авторов показали (Tiao and Xu. Robustness of maximum likelihood estimates for multi-step predictions: The exponential smoothing case, Xia and Тong. Feature matchingin time series modeling), что использование критерия \eqref{eq:estim_MSE1tohsteps} в условиях, когда прогнозная модель специфицирована неверно, позволяет увеличить точность прогнозов. А, как мы знаем из параграфа «Системы и модели«, все модели неверны. Даже в том случае, если модель по какой-то причине верна (спасибо божественному вмешательству), могут возникнуть проблемы из-за конечности выборки (у нас же не может быть бесконечного числа наблюдений, так ведь?), по которой оценивается модель. В этом случае использование методов оптимизации на основе многошаговых прогнозов так же позволяет увеличить точность прогнозов.
В качестве небольшого заключения заметим, что использование целевой функции \eqref{eq:estim_MSE1tohsteps} прекрасно вписывается в идею об инерционности прогнозируемых систем и сроках прогнозирования: мы должны учитывать, что разные объекты прогнозирования ведут себя по-разному, а значит и при построении прогнозных моделей для них нужно учитывать их инерционность и соотносить её со сроком. Легче всего это сделать, если использовать целевую функцию, соответствующую поставленной прогнозной задаче.
Траекторная функция правдоподобия
В этом подпараграфе мы поговорим об идеи, которая посетила меня после очередного обсуждения с Никосом на общую тему по прогнозированию. Мы уже упоминали ранее, что в целевой функции \eqref{eq:estim_MSE1tohsteps} почти всё хорошо, ей только не хватает статистического обоснования. Для того, чтобы получить такое обоснование, было решено обратиться к функции правдоподобия. Полученный метод мы назвали «Траекторной функцией правдоподобия» («Trace likelihood»). Слово «траекторный» здесь должно обозначать, что при использовании целевой функции мы делаем прогноз по некой траектории (от 1 до \(h\) шагов вперёд). На русский это перевести тяжело, но слово «траектория» звучит достаточно близко к тому, что имеется в виду.
Итак. Вместо того, чтобы смотреть на распределение одной ошибки (например, ошибки на одно наблюдение вперёд), мы решили изучить совместное распределение нескольких ошибок: от \(e_{t+1|t} \text{ до } e_{t+h|t}\). Для того, чтобы сделать это нам нужно обратиться уже не просто к дисперсии ошибок (как мы это делали в предыдущем параграфе), а к ковариационной матрице, которая содержит в себе как дисперсии соответствующих ошибок, так и ковариации между ними:
\begin{equation} \label{eq:estim_tracelikelihoodSigma}
{\Sigma} =
\begin{pmatrix}
\sigma_1^2 & \sigma_{1,2} & \dots & \sigma_{1,h} \\
\sigma_{1,2} & \sigma_2^2 & \dots & \sigma_{2,h} \\
\vdots & \vdots & \ddots & \vdots \\
\sigma_{1,h} & \sigma_{2,h} & \dots & \sigma_h^2
\end{pmatrix} ,
\end{equation}
По диагонали \eqref{eq:estim_tracelikelihoodSigma} как раз располагаются дисперсии, а вне диагонали — ковариации между ошибками. Заметим, дисперсии и ковариации, указанные в этой матрице, — на самом деле не простые, а условные. То есть они нам показывают, какой будет соответствующая величина при заданных параметрах модели. Условные ковариации здесь показывают, коррелируют ли друг с другом наши прогнозные значения. То есть, если прогноз \(t+1\) оказался выше фактического значения, окажется ли прогноз на шаг \(t+2\) так же в среднем выше? Фактически введя такую матрицу, мы отказались от популярного среди статистиков допущения об отсутствии автокорреляции остатков в модели (об этом мы поговорим в следующих главах).
Имея ковариационную матрицу \eqref{eq:estim_tracelikelihoodSigma}, мы теперь можем рассмотреть совместное распределение ошибок.
Сделаем небольшое лирическое отступление. Все знают, как выглядит функция распределения одной величины, но не все понимают, что такое совместное распределение. Обратимся для примера к нашему любимому нормальному распределению. Мы уже кратко обсуждали функцию нормального распределения в «Статистическом анализе«, вот её изображение для напоминания:
Нормальное распределение
Функция нормального распределения — это своеобразный двумерный колокол, который показывает, с какой вероятностью могут появиться те или иные значения (лежащие в заданных промежутках). Так вот, когда мы говорим о многомерной функции нормального распределения, то мы просто этот плоский колокол трансформируем в многомерный. Например, для двух переменных нужно представить себе обычный трёхмерный колокол. Если он ровный, то он будет характеризовать такое распределение, в котором две случайные величины независимы. График именно такой зависимости (бессовестно стащенный мною из Википедии) и показан на следующем рисунке:
Многомерное нормальное распределение. Источник картинки — википедия.
Если же переменные связаны друг с другом, то этот колокол будет деформирован и вытянут по какой-нибудь диагонали. Понимая как выглядит функция распределения двух переменных, вы теперь легко можете себе представить функцию совместного нормального распределения трёх, четырёх и любого другого не отрицательного количества переменных.
Итак, по аналогии с функцией правдоподобия для одной переменной рассмотрим вероятность появления сразу нескольких фактических значениях при заданных параметрах:
\begin{equation} \label{eq:estim_traceprobability}
P(y_{t+1 |t}, y_{t+2 |t}, \cdots, y_{t+h |t} | \theta, \Sigma) ,
\end{equation}
где \(y_{t+1 |t}\) — условное фактическое значение на наблюдении \(t+1\) при наличии информации на наблюдении \(t\). Для простоты ряд значений \( y_{t+1 |t}, y_{t+2 |t}, \cdots, y_{t+h |t}\) можно объединить в вектор, который мы будем обозначать через \(Y_t\). Тогда на основе \eqref{eq:estim_traceprobability} правдоподобие может быть записано так:
\begin{equation} \label{eq:estim_traceprobabilitylikelihood}
L(\theta, \Sigma | Y_t) = P(Y_t | \theta, \Sigma) .
\end{equation}
Ну, а теперь рассмотрим вероятность получения таких значений по всей выборке из \(T\) наблюдений:
\begin{equation} \label{eq:estim_generaltracelikelihood}
L(\theta, \Sigma | Y) = \prod_{t=1}^T P(Y_t | \theta, \Sigma) ,
\end{equation}
где \(Y = \begin{pmatrix}Y_1 \\ Y_2 \\ \dots \\ Y_T \end{pmatrix}\).
В принципе, уже можно максимизировать функцию правдоподобия \eqref{eq:estim_generaltracelikelihood}, только, опять же, как мы это любим делать, сделаем предположение о распределении наших случайных величин. Пусть оно будет многомерным нормальным (как ранее и обсуждалось):
\begin{equation} \label{eq:estim_generaltracelikelihood1}
L(\theta, \Sigma | Y) = \prod_{t=1}^T \left[ \left( 2 \pi \right) ^ {-\frac{h}{2}} | \Sigma |^{-\frac{1}{2}} \exp \left( -\frac{1}{2} {e_t}’ \Sigma^{-1} {e_t} \right) \right] .
\end{equation}
Однако, здесь надо иметь в виду, что \( {e_t} = \begin{pmatrix} \epsilon_{t+1} \\ \epsilon_{t+2} \\ \vdots \\ \epsilon_{t+h}\end{pmatrix}\) — это вектор условных ошибок, а не единичное значение.
Линеаризуем функцию \eqref{eq:estim_generaltracelikelihood1} с помощью натурального логарифма:
\begin{equation} \label{eq:estim_logtracelikelihood}
l(\theta, \Sigma | Y) = -\frac{T}{2} \left( h \ln(2 \pi) + \ln | \Sigma | \right) -\frac{1}{2} \sum_{t=1}^T \left( {e_t}’ \Sigma^{-1} {e_t} \right)
\end{equation}
Теперь по аналогии с тем, как мы оценивали дисперсию на основе функции правдоподобии ранее, мы можем оценить ковариационную матрицу \eqref{eq:estim_tracelikelihoodSigma}:
\begin{equation} \label{eq:estim_tracelikelihoodSigmaest}
\hat{\Sigma} = \frac{1}{T} \sum_{t=1}^T {e_t} {e_t’} .
\end{equation}
Вставляя \eqref{eq:estim_tracelikelihoodSigmaest} в \eqref{eq:estim_logtracelikelihood}, после ряда элементарных преобразований, получим логарифмированную концентрированную функцию правдоподобия:
\begin{equation} \label{eq:estim_logtracelikelihoodconcentrated}
l(\theta | Y) = -\frac{T}{2} \left( h \ln(2 \pi e) + \ln |\hat{\Sigma}| \right)
\end{equation}
Её теперь можно использовать при выборе прогнозных моделей (на основе информационных критериев), а так же для расчёта ковариаций и дисперсий коэффициентов. Но нас пока всё ещё интересует максимизация этого критерия. Окинув \eqref{eq:estim_logtracelikelihoodconcentrated} взглядом, мы придём к выводу о том, что максимум функции правдоподобия будет соответствовать минимуму следующей функции:
\begin{equation} \label{eq:estim_GV}
TFL = |\hat{\Sigma}| .
\end{equation}
TFL расшифровывется как «Trace Forecast Likelihood» (Функция правдоподобия траекторных прогнозов). Теперь самое сложное — понять, что же значит минимизация определителя ковариационной матрицы. Для этого мы вспомним формулу коэффициента корреляции, упомянутую ранее в предыдущих параграфах:
\begin{equation} \label{eq:estim_correlcoefforGV}
r_{x,y} = \frac{\sigma_{x,y}}{\sigma_x \sigma_y} .
\end{equation}
Из формулы \eqref{eq:estim_correlcoefforGV} можно выразить ковариации:
\begin{equation} \label{eq:estim_covarsforGV}
\sigma_{x,y} = r_{x,y} \sigma_x \sigma_y .
\end{equation}
Эти ковариации ничто не мешает подставить в нашу ковариационную матрицу \eqref{eq:estim_tracelikelihoodSigma} следующим образом:
\begin{equation} \label{eq:estim_tracelikelihoodSigmacorrel}
\Sigma =
\begin{pmatrix}
\sigma_1^2 & r_{1,2} \sigma_1 \sigma_2 & \dots & r_{1,h} \sigma_1 \sigma_h \\
r_{1,2} \sigma_1 \sigma_2 & \sigma_2^2 & \dots & r_{2,h} \sigma_2 \sigma_h \\
\vdots & \vdots & \ddots & \vdots \\
r_{1,h} \sigma_1 \sigma_h & r_{2,h} \sigma_2 \sigma_h & \dots & \sigma_h^2
\end{pmatrix} .
\end{equation}
Эта же прекрасная квадратная матрица \eqref{eq:estim_tracelikelihoodSigmacorrel} может быть представлена как произведение следующих трёх матриц:
\begin{equation} \label{eq:estim_tracelikelihoodSigmacorrel2}
\Sigma =
\begin{pmatrix}
\sigma_1 & 0 & \dots & 0 \\
0 & \sigma_2 & \dots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \dots & \sigma_h
\end{pmatrix}
\begin{pmatrix}
1 & r_{1,2} & \dots & r_{1,h} \\
r_{1,2} & 1 & \dots & r_{2,h} \\
\vdots & \vdots & \ddots & \vdots \\
r_{1,h} & r_{2,h} & \dots & 1
\end{pmatrix}
\begin{pmatrix}
\sigma_1 & 0 & \dots & 0 \\
0 & \sigma_2 & \dots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \dots & \sigma_h
\end{pmatrix}
\end{equation}
По середине этого сандвича \eqref{eq:estim_tracelikelihoodSigmacorrel2} теперь находится корреляционная матрица, а по бокам — простые диагональные матрицы со стандартными отклонениями. Если теперь взять определитель этого бутерброда, то (следуя одному из базовых свойств определителя) получится следующее:
\begin{equation} \label{eq:estim_tracelikelihoodSigmacorrel3}
|\Sigma| =
\prod_{j=1}^h \sigma_j^2
| R |
\end{equation}
где \(R = \begin{pmatrix}
1 & r_{1,2} & \dots & r_{1,h} \\
r_{1,2} & 1 & \dots & r_{2,h} \\
\vdots & \vdots & \ddots & \vdots \\
r_{1,h} & r_{2,h} & \dots & 1
\end{pmatrix}\) — упомянутая выше корреляционная матрица.
Логарифм этого безобразия нам даст чуть более приятное безобразие:
\begin{equation} \label{eq:estim_tracelikelihoodSigmacorrel4}
\ln |\Sigma| =
\sum_{j=1}^h \ln \sigma_j^2 +
\ln | R |,
\end{equation}
Анализируя его, можно заметить, что при минимизации обобщённой дисперсии \eqref{eq:estim_GV} мы фактически минимизируем логарифмы дисперсий и одновременно с этим уменьшаем логарифм определителя корреляционной матрицы. Если с первым всё более-менее понятно, то со вторым — не совсем. Для того, чтобы понять, что это значит, рассмотрим один простой пример. Предположим, что нас интересует прогноз на срок \(h=2\). Корреляционная матрица в \eqref{eq:estim_tracelikelihoodSigmacorrel4} в этом случае становится размером 2 на 2. Из курса высшей математики мы знаем, что определитель такой матрицы будет равен:
\begin{equation} \label{eq:estim_logtracelikelihoodconcentratedCF2}
| R | = 1 \cdot 1 — r_{1,2} \cdot r_{1,2} = 1 — r_{1,2}^2
\end{equation}
Минимум этого определителя будет соответствовать ситуации, когда корреляция \(r_{1,2} = 1\) или \(r_{1,2} = -1\). То есть минимизируя определитель, мы увеличиваем существующую корреляцию между прогнозными ошибками. Здесь стоит оговориться: если между ними корреляции нет, то и увеличивать особо нечего. То есть этот эффект будет наблюдаться только в ситуациях, в которых уже есть проблема коррелированности прогнозных ошибок.
Обобщая этот результат и обращаясь к курсу многомерного статистического анализа, можно заключить: использование критерия \eqref{eq:estim_GV} будет приводить к увеличению корреляции между теми ошибками, которые уже итак коррелируют друг с другом. Не самое приятное свойство, надо сказать. Причём, увеличив корреляцию, мы добавляем проблем к своей модели, так как эту корреляцию нужно как-то учесть. В общем, пока, как быть с этой ситуацией, не вполне понятно.
Однако у любой проблемы есть решение. И самое простое решение в нашем случае — это проигнорировать проблему. Давайте предположим, что корреляции в нашей матрице нулевые, то есть матрица \( R \) диагональная (и к тому же единичная), и имеет следующий вид:
\begin{equation} \label{eq:estim_tracelikelihoodSigmaTLV}
R = I_h =
\begin{pmatrix}
1 & 0 & \dots & 0 \\
0 & 1 & \dots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \dots & 1
\end{pmatrix} ,
\end{equation}
Её определитель будет равен единице, а значит критерий \eqref{eq:estim_tracelikelihoodSigmacorrel4} трансформируется в нечто более простое и понятное:
\begin{equation} \label{eq:estim_TLV}
GMSTFE = \sum_{j=1}^h \ln \sigma_j^2 ,
\end{equation}
GMSTFE расшифровывается как «Geometric Squared Trace Forecast Error» (Геометрическая квадратная траекторная ошибка прогноза). Название, возможно, не самое красивое, но отображает примерный смысл происходящего.
В чём же преимущества в использовании целевых функций \eqref{eq:estim_GV} и \eqref{eq:estim_TLV} по сравнению, например, с функцией \eqref{eq:estim_MSE1tohsteps}? Во-первых, они статистически обоснованы, что позволяет на основе них рассчитывать значение функции правдоподобия, а значит и легко использовать информационные критерии при выборе моделей. Во-вторых, эти функции дают информацию о распределении коэффициентов (упомянутые всуе уже несколько раз дисперсии коэффициентов). Ну, и в-третьих, у них есть полезное свойство: логарифм позволяет изменить шкалу величин.
Мы уже упоминали ранее, что дисперсия одношаговых ошибок обычно меньше дисперсии ошибок на несколько шагов вперёд. Логарифмирование делает их близкими друг к другу. Пусть, например, у нас есть дисперсии ошибок на 1 шаг и 6 шагов вперёд: \(\sigma_1^2 = 500 \), \(\sigma_6^2 = 3000 \). Возьмём их логарифмы: \(\ln (\sigma_1^2) \approx 2.7 \), \(\ln(\sigma_6^2) = 3.5 \). Видно, что величины стали ближе к друг другу, хотя изначально разница между ними составляла 6 раз. Это означает, что при минимизации \eqref{eq:estim_GV} или \eqref{eq:estim_TLV} не будет эффекта перекрывания дисперсий, который наблюдался в \eqref{eq:estim_MSE1tohsteps} — дисперсии краткосрочных ошибок будут минимизироваться наравне с дисперсиями долгосрочных.
Ну, и в качестве заключения этого параграфа пара мыслей.
Первая: использование методов \eqref{eq:estim_GV} и \eqref{eq:estim_TLV} имеет смысл лишь в моделях авторегрессий и экспоненциального сглаживания.
Вторая: обе функции прекрасно вписываются в идею с инерционностью систем и сроками прогнозирования. Если вы собираетесь сделать долгосрочный прогноз продаж шляп в Петербурге, то использовать критерий, учитывающий исключительно краткосрочные ошибки, неразумно. Лучше обратиться к \eqref{eq:estim_TLV}, \eqref{eq:estim_GV}, ну или на худой конец к \eqref{eq:estim_MSE1tohsteps} или \eqref{eq:estim_MSEhsteps}.