CC BY
25
Методика Интеграции Формальных Методов Прогнозирования Временных Рядов в Data
Assimilation
Ю.С. Тимошенкова, С.В. Поршнев, Н.Т. Сафиуллин
Аннотация—В статье описан разработанный авторами способ интеграции в метод Data Assimilation (DA) формальных методов прогнозирования временных рядов (ВР) (авторегресии-проинтегрированного скользящего среднего, сингулярного спектрального анализа, метода группового учета аргумента, нейронных сетей с долговременной памятью), применяемых в тех случаях, когда неизвестна математическая модель динамической системы, породившей данный ВР (например, ВР, составленные из значений экономических показателей). Работоспособность предложенного способа интеграции проиллюстрирована на примере прогнозирования ВР «Air Passengers» методом DA на основе ансамблевого фильтра Калмана, в который интегрирован метод ARIMA. Получены оценки точности прогнозов значений выбранного ВР, вычисленные с помощью ARIMA и с помощью предложенного метода. Обсуждаются преимущества и недостатки предложенного способа интеграции метода DA и формальных методов прогнозирования ВР, а также направления его дальнейшего совершенствования.
Ключевые слова—Анализ временных рядов, прогнозирование временных рядов, Data Assimilation.
I. Введение
Задача прогнозирования одномерного временного ряда (ВР), под которым понимают значения некоторого процесса, порожденного изучаемой динамической системой (ДС), зафиксированные в упорядоченные в порядке возрастания моменты времени, является актуальной в различных областях человеческой деятельности, например, в экономике [1], метеорологии, климатологии [2], океанографии [3] и др. В общем случае данная задача состоит в предсказании на основе известных значений анализируемого ВР его значений в последующие моменты времени.
Множество задач прогнозирования ВР в зависимости от типа ВР можно разделить на следующие классы.
Статья получена 17 января 2022
Тимошенкова Юлия Сергеевна, Уральский Федеральный Университет, Институт радиоэлектроники и информационных технологий-РТФ, (e-mail: julia.timoshenkova@urfu.ru).
Поршнев Сергей Владимирович, Уральский Федеральный Университет, Институт радиоэлектроники и информационных технологий-РТФ, (e-mail: s.v.porshnev@urfu.ru).
Сафиуллин Николай Тахирович, Уральский Федеральный Университет, Институт радиоэлектроники и информационных технологий-РТФ, (e-mail: n.t.safiullin@urfu.ru).
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 20-37-90006
1. Прогнозирование ВР, порожденных ДС, для
описания эволюции которых существуют те или иные математические модели.
2. Прогнозирование ВР, порожденных ДС, не
имеющих математических моделей вследствие их сложности.
Задачи первого класса возникают в автоматизированных системах управления различными техническими системами, например, ракетно-космической техники [4], задачи второго класса - в экономике, например, при прогнозировании курсов валют [5].
Соответственно, изучают два класса методов прогнозирования ВР:
1) методы прогнозирования ВР, опирающиеся на строгое математическое описание ДС, породившей данный ВР (неформальные методы);
2) методы прогнозирования ВР, опирающиеся только на сам ряд исходных наблюдений, без учета наличия или отсутствия информации об исходной ДС, породившей данный ВР (формальные методы);
К неформальным методам прогнозирования ВР относится, например, метод фильтра Калмана [6] (во всех его известных модификациях), к формальным методам прогнозирования ВР - метод авторегрессии и скользящего среднего (Autoregressive Integrated Moving Average, ARIMA) [7]; метод сингулярного спектрального анализа (Singular spectrum analysis, SSA), в русскоязычной научной литературе известный как метод «Гусеница» [8]; метод группового учета аргумента (МГУА), в англоязычной литературе - Group Method of Data Handling (GMDH) [9]; метод прогнозирования на основе использования рекуррентных нейронных сетей (Recurrent Neural Network, RNN) [7].
На практике решается, в том числе, задача прогнозирования таких ВР, у которых известны как спрогнозированные, так и соответствующие фактические (измеряемые) значения (например, краткосрочное прогнозирование погоды). Наличие точных и спрогнозированных значений ВР, синхронизированных во времени, обеспечивает возможность их сравнения и при необходимости последующей корректировки используемой модели ВР и/или уточнения спрогнозированных значений данного ВР в последующие моменты времени. Описанный подход реализован, например, в методе ассимиляции данных (Data Assimilation, DA) [4], относящийся к
классу неформальных методов прогнозирования ВР.
Напомним, что сегодня известны несколько модификаций БЛ, в том числе: 1) модификация, основанная на вычислении прогнозируемых значений ВР через построение функционала, обеспечивающего близость значений ВР, рассчитанных на основе математической модели динамической системы, породившей данный ВР, и наблюдаемых значений ВР; 2) модификация, основанная на использовании статической теории оценивания и фильтра Калмана, в которой ВР рассматривается как аддитивная смесь информативного сигнала и случайного процесса (шума) с известными статистическими свойствами.
Примеры успешного использования БЛ для прогнозирования ВР приведены в [10]—[13]. Однако, как показал анализ многочисленных работ, посвященных формальным методам прогнозирования ВР (см., например, [13] и др.), данные методы совместно с БЛ не использовались. В этой связи разработка способов интеграции формальных методов прогнозирования ВР и БЛ является актуальной задачей.
В статье описан разработанный авторами способ интеграции формальных методов прогнозирования ВР в метод БЛ, основанный на использовании статистической теории оценивания и фильтра Калмана, приведены примеры использования для прогнозирования реальных ВР и получены оценки точности прогноза.
II. Материалы и Методы Работы Л. Классический БЛ
В методе БЛ полагается, что эволюция исследуемой ДС, характеризующейся в дискретном времени (к
вектором состояния Хк, описывается следующей системой уравнений [2], [12]:
x;+i = м (x; > fk ) + w = H ( xk > fk )+
(1)
Изначально метод БЛ, обсуждаемый далее, был разработан для линейной ДС, в математической модели которой операторы Ми Н в (1) являются матрицами, поэтому
м (х л ) = м ■ X (Ч ) = м ■ Х = мк, н (Хк Л ) = Н ■ Хк (Гк ) = Н ■ Хк = Нк, а , Ук - выборки из нормального белого шума с известными значениями ковариационных матриц ошибок Qк, Як.
На этапе прогнозирования на основе известного в момент времени tк вектора состояния ДС хк и вектора наблюдений ук вычисляют скорректированное значение вектора состояния ДС по формуле:
xa=xf+k:( ёк ),
(2)
где
ек = ук - Н (х[ ) - отклонение реального вектора
наблюдений ук от спрогнозированного вектора наблюдений, полученного на основе вектора состояния ДС Хк , в линейном случае: ~ёк = ук - Н ■ х( = ук - Нк;
К * - передаточный коэффициент фильтра Калмана, вычисляемый по формуле:
К; = Pkl (HkPkl + Rk )-
(3)
Здесь Pf - ковариационная матрица ошибок прогноза, вычисляемая по формуле:
Pf = мк (I - к;., нк-, ) F/-,mT + Qk
(4)
здесь I - единичная матрица.
Также на данном этапе вычисляют ковариационную матрицу ошибок прогноза Р£ :
pa=( I - к* нк ) ppf.
(5)
где Хк, к = 0,1,... - вектор состояния динамической системы в моменты времени /к, к = 0,1,...; М -оператор перехода (математическая модель, описывающая динамику системы); Хк+1 -спрогнозированное значение вектора состояния ДС в момент времени (к+х, ук - вектор наблюдений,
Н (Хк, tk) - оператор наблюдения, связывающий вектор состояния ДС с вектором наблюдений; м>к, Ук - ошибки
модели и наблюдения, соответственно.
Первое уравнение в (1) называется уравнением прогноза, второе уравнение - уравнением наблюдений. Соответственно, алгоритм БЛ разделяется на два этапа: этап прогнозирования и этап корректировки прогноза. Отметим, что корректировка вектора Хк в методе БЛ может осуществляться не на каждом шаге прогноза, но в моменты времени tk+т, tk+2т, ..., т - целое число, равное или большее единицы.
На этапе корректировки прогноза с учетом вектора наблюдений ук в^гчисляют уточненное значение
вектора состояния ДС Хк по формуле:
ха = Х{ + К; (Ук - Н (Х{)) = Х{ + к; (у к - н к), (6)
а также новую ковариационную матрицу ошибок прогноза Р/+1:
Pf = M P aMT + Q к+1 1v1 к+1гк 1v1 к+1 т \ск ■
(7)
используемую далее для вычисления xkf+1 по формуле:
7f -
=м ( xp ).
(8)
Затем повторяется описанная выше последовательность действий.
Блок-схема алгоритма БЛ, описанного, следуя [2], [14], представлена на рисунке 1.
Поясним смысл оператора модели и оператора наблюдения в (1) на примере линейного гармонического осциллятора (например, точечного груза массой т, прикрепленного к одному из концов пружины с жесткостью к, второй конец которой закреплен неподвижно), совершающего свободные колебания, который описывается дифференциальным уравнением (ДУ) второго порядка:
r +ю r = 0,
(9)
где г - смещение точечного груза относительно равновесного положения, ю = ^к/т - циклическая частота колебаний гармонического осциллятора, с начальными условиями г (0) = г0, Г (0) = у0 ; и на примере системы дифференциальных уравнений первого порядка (СДУ), предложенной Э.Н. Лоренцом для описания процесса конвекции в плоском двумерном слое морской воды [15], [16]:
X = ст(y - x) y = x(r - z) - y, z = xy - bz
(10)
z1 z2
_ z2 _ -ю2 z1
[ Z1
Используем в (11) конечно-разностную аппроксиммацию производных первого порядка точности:
С
11ачало
Задать параметры метода. DА
)
Г
Цикл j от len до len + m
Выполнить прогноз точки к с использованием математической модели (1)
т
В точке к получить данные наблюдений системы ук
Сравнить прогноз и наблюдение и точке к согласно (2)
Скорректировать прогнозируемое значение согласно (4) (Ci)
V
Обновись матрицы ковариаций согласно результатам коррекции (7)
I
Увеличитьj на [
С
поразить полученные результаты про г
Конец
пноза j
)
где ст = у/к - число Прандтля, г = Ка/Кас -
нормированное число Рэлея, Ь = 4^ (1 + а2) -
геометрический фактор, х - интенсивность конвекции, у - разность температур между восходящим и нисходящим потоками; г - отклонение вертикального температурного профиля от линейного. СДУ (10) также возникает в следующих математических моделях: конвекции в замкнутой петле; вращения водяного колеса; одномодового лазера; диссипативного гармонического осциллятора с инерционной нелинейностью. Напомним, что СДУ (10) замечательна тем, что при ст = 10, г = 28, Ь = 8/3 ее решение, получившее название аттрактор Лоренца, оказывается хаотическим.
Для записи в явном виде оператора эволюции М гармонического осциллятора представим ДУ второго порядка (10) в виде системы ДУ первого порядка:
(11)
где = г, = ¿1, г2 = г1 = г =-юг = -юг1. Из (11) видно, что в пространстве состояний линейный гармонический осциллятор описывается вектором
Рис. 1 - Классический алгоритм метода БЛ: т - длина прогноза, к - номер отсчета прогноза, 1еп - длина вектора наблюдений
Z1 & )" Z2 (h )_
[ Z1 ]k+1 -[ Z1 ]k T
[ Z2 ]k+1 -[ z2 ]k
Г[ z, ] ]
L U k+1 =
[ z2 ]
X 2J k+1 _
" 1 t] Г[ zi 1]
= L U k
-Ю2 T 1 _[z2 ]k _
. (12)
M =
где т - шаг интегрирования СДУ(11). Тогда можно записать модель линейного гармонического осциллятора в дискретной форме:
[]к +[22 ]к т _-®2 Т[ ]к +[ 22 ]к _ Из (12) видно, что искомый оператор эволюции линейного грамонического осциллятора есть матрица
" 1 т" -ю2т 1_
В тех случаях, когда проводятся измерения только координаты или скорости линейного гармонического осциллятора, вектор наблюдений ук во втором уравнении системы (1) вырождается в один из скаляров:
Ук =[ ]к, У к =[ г2 ]к ,
соответственно. Следовательно, в рассматриваемых случаях оператор наблюдений Н представляет собой одну из вырожденных матриц:
" 1 0"
0 0
H =
H _
0 0 0 1
соответсвенно.
В тех случаях, когда проводятся одновременные измерения и координаты и скорости гармонического
1 * И"
осциллятора, вектор наблюдений будет yk _
[]
Соответственно, оператор наблюдений Н представляет собой единичную матрицу.
В тех случаях, когда измеряется ускорение гармонического осциллятора г ^) = ¿2 (t), вектор ук во втором уравнении системы (1) вырождается в скаляр ук = Н ([ 2Х \к, [ ]к, [ ]к). Принимая во внимание
связь между []к, []к, [¿2 ]к, понятно, что в
рассматриваем случае оператор наблюдений Н является оператором вычисления второй производной от 2Х (t), для вычисления которого можно использовать следующую конечно разностную аппроксимацию:
H ([L ) = [ -2 ]k -Т[-2 ]k-1.
выбор обусловлен тем, что, с одной стороны, он имеет относительно невысокую в сравнении с другими модификациями вычислительную сложность, а с другой стороны, обеспечивает приемлемую точность прогнозирования [12]. Блок-схема БЛ, в котором используется ЕпКР-фильтр, представлена на рисунке 2.
Из рисунка 2 видно, что в данной модификации метода БЛ в каждый момент времени tk генерируется ансамбль наблюдений:
—(ensemble) — , —
Ук ,, _ Уь + u>
. _о jj{ensemble)
(14)
где ц. - шум с известным законом распределения. Далее, используя случайные значения реализации —УететЫе), в соответствие с (6), вычисляется ансамбль
Ук.
скорректированных значений xt
i(eensemble^
а в качестве
скорректированного значения прогноза - усредненное
по ансамблю x.
—a^ensemble)
значение х, _ mean
—a^ensemble^ \
Хк j у
При этом вместо корреляционной матрицы ошибок наблюдения используется ее несмещенная оценка, вычисляемая по формуле:
Для нахождения в явном виде оператора эволюции м для СДУ (10), используем конечно разностную аппроксиммацию первых производных
Т
У,+1 - У,
Т
+1 -
где т - шаг интегрирования СДУ (10), запишем модель Лоренца в дискретной форме:
xt+1 1 -стт стт 0 X
У,+1 _ r Т 1 -Т -X, Т У, . (13)
3+1 _ Уt т 0 1 - b т 3 _
Из (13) видно, что при использовании классического метода БЛ оператор наблюдения в (1) представляет собой матрицу:
M _
1 -стт стх 0 r Т 1 - Т -X, y, т 0 1 - b т
Оператор наблюдений H при предположении о том, что одновременно наблюдаются значения всех трех переменных x, y, z, - оказывается единичной матрицей.
Отметим, что сегодня разработано множество различных модификаций фильтра Калмана, используемых в DA, в том числе: EKF (Extended Kalman filter), EnKF (Ensemble Kalman Filter), UKF (Unscented Kalman Filter) и др. [2], из которых далее в работе мы используем EnKF (ансамблевый фильтр Калмана). Его
Рис. 2 - Схема алгоритма метода БЛ с использованием ЕпКР-фильтра: т - длина прогноза, к - номер отсчета прогноза, 1еп - длина вектора наблюдений, N - размер ансамбля фильтра
R _
u0u0 U 00U\ ' U0UM
UjUj u1um
Здесь чертой сверху обозначена процедура усреднения по ансамблю реализаций ц, а вместо произведений
Р/нТ , НкРк 'Нтк - используются их статистические оценки:
PfHT _ ■
гк 11 к
- H
x - mean
( Xlit ))( H ( Xfi, )- mean ( у{ , ))T
H PfHT _ 11krk 11 к
Pf _
Z(xf-m -Ук)(xf-Ук)T,m _ 01. (15)
_ F ( xk ,ak> ,k ) + i
У к _ X + vk,
к-( p-1)J
= T+J Z( H (Xf,,)-mean (% ))•
• ( xf - mean (xI)) ( н ( Xfki) - mean ( у I ))r ,
и вместо ковариационной матрицы ошибок прогноза -ее оценка, рассчитываемая по формуле:
где хк =(хк, Хк -1,хк_р )к , ак = (ак, а к-¡, Р - функция, определяемая выбранным формальным методом прогнозирования ВР; ук - наблюдаемое значение вектора состояния ДС, здесь, оператор наблюдения, связывающий Хк с наблюдаемым значением ук, Н = 1, , Ук - ошибки модели и
наблюдения, соответственно. Вектор параметров математической модели (16) при использовании формальной модели в методе БЛ идентифицируется на основе использования к, к -1,..., к - (р +1) значений ВР.
Блок-схема алгоритма интеграции формальных методов прогнозирования ВР в метод БЛ (17), основанном на использовании ЕпКР-фильтра Калмана, представлен на рисунке 3.
Из рисунка 3 видно, что данный алгоритм (далее алгоритм интеграции) реализуется выполнением следующей последовательности действий:
1. Выбор формального метода прогнозирования.
2. Задать параметры метода БЛ и формального метода прогноза.
Результаты сравнительного анализа качества прогнозирования классическим методом БЛ и методом БЛ, использующим ЕпКР-фильтр, представлены в [17].
В. Интеграция формальных методов
прогнозирования ВР в метод БЛ
При решении задачи прогнозирования ВР, порожденных ДС, для которых на данный момент не существует соответствующих математических моделей, с помощью формальных методов прогнозирования, полагают, что
хк+1 f ( хк, xk-1 '...' xk - p, а1 '...' ap+1)
(16)
где Хк+1 - прогнозируемое значение; Р - функция, определяемая выбранным формальным методом прогнозирования ВР; Хк,Хк-1,...,Хк-р - предыдущие
значения ВР, а1,...,ар+1 - параметры модели.
Сравнивая (16) и уравнение прогноза в (1) видим, что с точностью до аддитивного члена (ошибки
прогноза), данные уравнения имеют одинаковую структуру, поэтому функция
по сути, аналогична
(1). Следовательно,
использование (16) эквивалентно переходу от описания ДС в пространстве состояний к ее описанию в р-мерном пространстве координатном пространстве. Используя отмеченную аналогию, можно сразу записать систему уравнений, объединяющую метод БЛ и формальные методы прогнозирования ВР в виде:
F (xk, Xk-1,..., xk-p ,a\,...,am ) ,
оператору M ( xk, ,к )
Рис. 3 - Алгоритм интеграции формальных методов прогнозирования ВР в метод БЛ
3. Получить ансамбль прогнозов в точке с помощью формального метода (в рамках работы — это прогноз с использованием АР) и добавления шума.
4. Получить данные наблюдений системы в прогнозируемый момент времени.
5. Рассчитать матрицы ковариации ошибок прогноза.
6. Выполнить коррекцию прогноза путем применения матриц ковариации и передаточного коэффициента фильтра Калмана K для каждой точки ансамбля.
7. Обновить каждый член ансамбля прогноза и матрицы ковариаций в соответствии с полученным значением K.
8. Выполнить пункты 3-7 для каждой прогнозируемой точки.
9. Оценить точность полученного прогноза.
В данном решении параметры используемых моделей задаются один раз в начале и не меняются со временем выполнения прогнозирования.
C. Описание прогнозируемого ВР
Тестирование разработанной методики проводилось на наборе данных «Air Passengers» [18], одного из популярных ВР, используемых для тестирования различных методов прогнозирования ВР [19]. Данный ВР составлен из ежемесячных объемов международных пассажироперевозок авиалиний США в период с января 1949 по декабрь 1960 гг., измеряемых в тысячах человек. Число членов ВР (длина ВР) - N = 144. Визуализация ВР в виде зависимости «объема ежемесячных пассажироперевозок от времени» представлена на рисунке 4.
Анализ ВР, представленного на рисунке 4, позволяет высказать предположение о том, что он представляет собой линейную комбинацию отсчетов монотонно возрастающей (тренд) и периодической функций времени. Можно показать [20], что наилучшей аппроксимацией обсуждаемого тренда является линейная функция (см. Рисунок 5).
Разность между обсуждаемым ВР и линейным трендом представлена на рисунке 6.
Из рисунка 6 видно, что в данном ВР присутствует периодическая(с периодом 12 месяцев) и сезонная компоненты, ответственная за увеличение амплитуды
^HiOCO-ir-^T-H^tCOoi ОО^НООт-^ОО-—
месяц, год
Рис. 4 - Визуализация тренда и ВР «Air Passengers»: 1 - отсчеты ВР, 2 - тренд
^н Ю О СТ ^н Tf ОО OJ
ОО^ООт-нОО^!
месяц, год
Рис. 5 - Визуализация ВР «Air Passengers» с удаленным линейным трендом
периодической составляющей в период с июня по сентябрь, которые в [20] выделены на основе количественного анализа ВР «Air Passengers». Здесь наличие сезонной компоненты обусловлена тем, что в США большинство людей используют отпуска в летние месяцы. Также отметить, что анализируемый ВР ряд не имеет явных аномалий и выбросов, в связи с чем не требуется его предобработка с целью удаления выбросов.
В [21] обосновано, что данный временной ряд может быть описан одной из известных формальных авторегрессионных моделей, из которых, без потери общности полученных результатов, в нашем
исследовании будут использоваться AR-модели вида:
p
xt+1 = c + Хал+sk , (18) i =1
где c - константа, ек - белый шум,
ах,а2,...,ар - параметры модели (коэффициенты
авторегрессии). Напомним, что параметры AR-модели (18) находятся как решения уравнения Юла-Уокера [1], которое при определенных условиях сводится к методу наименьших квадратов.
III. Методика исследования и анализ результатов прогнозирования ВР «Air Passengers»
исследование было проведено в соответствие с методикой, реализующейся выполнением следующей
Ci о rH со ■4« ю 1- ОО Ci
■■tf ю ю ю ю ю 1Л ю
Ci Ci Ci Ci Ci Ci Ci Ci Ci
I—1 I—< !-1 1-1 Т—1 -—I I—1 1—1 T-t
H lO о Ы 1—1 Tji со
о о I-H о о ,—1 о о
месяц, год
Рис. 6 - Временной ряд «Air Passengers»
последовательности действии:
1. Разбиение прогнозируемого ВР на известные и
прогнозируемые наблюдения, число которых составили = 84 и Щ = 60 отсчетов,
соответственно.
2. Выбор порядка АЯ-модели р.
3. Идентификация параметров АЯ--модели (18) [1].
4. Прогнозирование 60 значений ВР (в период с
января 1956 по декабрь 1960 гг.) с помощью формальной АЯ-модели (18) и предложенной интегрированной методики.
5. Оценка точности прогноза на основе сравнения
спрогнозированных значений ВР и соответствующих реальных известных наблюдений с помощью критерия:
RMSEp = —
p Nf
^ (Yi+Ntest t+Ntes
1/2
(19)
где Щу - число прогнозируемых наблюдений,
У - спрогнозированные значения ВР, У - известные отсчеты ВР в прогнозируемых точках.
Были вычислены оценки значений ЯМБЕ для АЯ-моделей, порядок которых изменялся от 1 до 10. Для получения статистически значимых результатов размер ансамбля фильтра в ЕпКЕ выбирался равным 50. Далее с помощью АЯ-модели (18) и разработанного алгоритма интеграции были вычислены в соответствие с (17) зависимости дисперсии разностей между точными и спрогнозированными ЯМ5Е12 = / (р), представленные на рисунке 7.
Из рисунка 7 видно, что функция ЯМБЕг = / (р) при увеличении порядка АЯ-модели от 1 до 10 монотонно убывает от 94,4 до 43,2, а ЯМБЕ2 = /(р) оказывается возрастающей от 3,4 до 27,1. Диапазоны изменения зависимостей ЯМБЕ1 = / (р) для других порядков
АЯ-моделей представлены в таблице 2.
Для пояснения данного результата рассмотрим зависимости спрогнозированных значений от времени для случаев р = 1,10, представленные на рисунках 8,9.
Из рисунков 8, 9 видно, что при увеличении порядка АЯ-модели улучшается качество прогнозирования тренда, который оказывает определяющее влияние на
100
со £
СЧ
50
1
-■
К ~
1 2 3 4 5
Порядок АР
Рис. 7 - Зависимость ЯМБЕ от порядка АЯ-модели, используемой для прогнозирования: 1 - АЯ--модель, 2 - алгоритм интеграции
Таблица 1. Диапазоны изменения параметра RMSE
Порядок АР RMSE
Алгоритм интеграции ^R-модель
min max min max
1 2,75 3,381 92,35 94,40
2 15,32 17,24 89,43 91,74
3 13,33 14,98 69,41 72,87
4 13,00 14,69 65,15 67,83
5 18,14 21,92 57,01 59,40
6 22,94 23,07 57,27 58,60
7 20,74 24,51 55,91 57,86
8 24,51 26,80 42,61 45,36
9 24,97 26,96 43,012 45,63
10 26,51 27,11 40,66 43,17
10 20 30 40 Номер точки прогноза
50
60
Рис. 8 - Результат прогнозирования ВР «Air Passengers»
на основе AR-модели 1-го порядка: 1- истинные значения ВР (серая линия), 2 - AR--модель, 3 - алгоритм интеграции (пунктирная линия)
600
400
300
.__-î Г y ; \_
п / J* 1 \ t. n 1
•7 Y" г Х-^/ч,/ AT \ :\ J * Г ч V
10 20 30 40 Номер точки прогноза
50
60
Рис. 9 - Результат прогнозирования ВР «Air Passengers» на основе AR-модели 10-го порядка: 1 - истинные значения ВР, 2 - AR-модель, 3 - алгоритм интеграции
значения анализируемого ВР. Отмеченное свойство AR-модели объясняет монотонное убывание функции RMSE1 = f (p). Отметим, что AR-модель 10-го порядка
пытается отследить наличие периодической и сезонной составляющих анализируемого ВР, однако, точность ее прогнозирования оказывается недостаточно высокой. Здесь можно ожидать, что увеличить точность прогнозирования анализируемого ВР удастся при использовании метода SARIMA, разработанного для прогнозирования ВР с сезонной составляющей (Seasonal ARIMA) [22].) Таким образом, для повышения точности прогнозирования ВР «Air Passengers» с помощью (17) целесообразно использовать AR-модели высокого порядка. Однако увеличение порядка AR-модели
приводит к увеличению размерности системы уравнений Юла-Уокера и, соответственно, к пропорциональному увеличению числа необходимых известных пред-прогнозных точек, что не всегда возможно на практике. При использовании предложенного алгоритма интеграции, наоборот, наименьшая погрешность прогнозирования оказывается в случае использования А^-модели первого порядка. Данный результат, с нашей точки зрения, обусловлен тем, что для модели первого порядка уравнение прогноза и уравнение коррекции оказываются линейны относительно входящих в них переменных.
IV. Заключение
Предложен алгоритм интеграции формальных методов прогнозирования ВР в метод DA, в котором используется EnKF-фильтр Калмана, работоспособность которого подтверждена результатами прогнозирования ВР «Air Passengers». Проведен сравнительный анализ точности прогнозов, полученных с помощью формальной модели ВР, и алгоритма интеграции, результаты которого позволяют сделать следующие выводы:
1. На примере одномерного ВР продемонстрирована возможность достоверного прогнозирования его значений с помощью метода DA, в котором вместо математической модели, описывающей закон изменения во времени состояния ДС, породившей изучаемый ВР, используется его формальная модель ВР.
2. Точность прогнозирования ВР «Air Passengers», оцениваемая дисперсий разностей между точными и спрогнозированными значениями изучаемого ВР, с помощью оригинальной методики, основанной на интеграции формальных методов прогнозирования ВР и метода DA, оказалась выше аналогичной величины в случае прогнозирования на основе AR-модели ВР любого порядка.
Подтверждение работоспособности предложенной методики прогнозирования состояния ДС в пространстве состояний больших размерностей является предметом последующих публикаций.
Библиография
[1] G. E. Box, G. M. Jenkins, G. C. Reinsel, и G. M. Ljung, Time series analysis: forecasting and control. John Wiley & Sons, 2015.
[2] L. Nerger и др., «SANGOMA: Stochastic Assimilation for the Next Generation Ocean Model Applications EU FP7 SPACE-2011-1 project 283580».
[3] M. Kanamitsu, «Description of the NMC global data assimilation and forecast system», Weather and forecasting, т. 4, вып. 3, сс. 335-342, 1989.
[4] А. В. Костров и А. М. Мардюк, «Об одной математической задаче теории движения ракет», Журнал вычислительной математики и математической физики, т. 34, вып. 10, сс. 1427-1443, 1994.
[5] L. Kilian и M. P. Taylor, «Why is it so difficult to beat the random walk forecast of exchange rates?», Journal of International Economics, т. 60, вып. 1, сс. 85-107, 2003.
[6] G. Welch и G. Bishop, «An introduction to the Kalman filter», 1995.
[7] G. P. Zhang, «Time series forecasting using a hybrid ARIMA and neural network model», Neurocomputing, т. 50, сс. 159-175, 2003.
[8] N. Golyandina, V. Nekrutkin, и A. A. Zhigljavsky, Analysis of time series structure: SSA and related techniques. CRC press, 2001.
[9] А. Г. Ивахненко, Долгосрочное прогнозирование и управление сложными системами. Техшка, 1975.
[10] Y. Timoshenkova, S. Porshnev и N. Safiullin, «On the possibility of correction of the forecasting of the Lorenz attractor dynamic characteristics using experimental data and data assimilation», в Journal of Physics: Conference Series, 2018, т. 1053, с. 012004.
[11] Y. Timoshenkova, S. Porshnev, N. Safiullin и O. Ponomareva, «On the Possibility of the Forecast Correction Lorenz Systeme Use Particle Filter», в 2018 International Conference on Applied Mathematics and Computational Science, ICAMCS. NET 2018, 2018, сс. 27-30.
[12] P. L. Houtekamer и H. L. Mitchell, «Data assimilation using an ensemble Kalman filter technique», Monthly Weather Review, т. 126, вып. 3, сс. 796-811, 1998.
[13] M. Ghil и P. Malanotte-Rizzoli, «Data assimilation in meteorology and oceanography», Advances in geophysics, т. 33, сс. 141-266, 1991.
[14] B. K. W. Lahoz и R. Menard, Data assimilation. Springer, 2010.
[15] E. N. Lorenz, «Deterministic nonperiodic flow», Journal of atmospheric sciences, т. 20, вып. 2, сс. 130-141, 1963.
[16] М. Табор, Хаос и интегрируемость в нелинейной динамике. 2001.
[17] M. Roth, G. Hendeby, C. Fritsche, и F. Gustafsson, «The Ensemble Kalman filter: a signal processing perspective», EURASIP Journal on Advances in Signal Processing, т. 2017, вып. 1, сс. 1-16, 2017.
[18] «Air Passengers». https://kaggle.com/abhishekmamidi/air-passengers (просмотрено дек. 01, 2021).
[19] M. Bocquet, «Introduction to the principles and methods of data assimilation in geosciences», Notes de cours, Ecole des Ponts ParisTech, 2014.
[20] Главные компоненты временных рядов: метод «Гусеница» / Под ред. Д.Л.Данилова, А.А.Жиглявского. 1997. Просмотрено: дек. 20, 2021. [Онлайн]. Доступно на: https://rusneb.ru/catalog/001980_000024_RU_FESSL_MAIN_1331818796 268482813/
[21] А. Н. Тырсин, «Построение моделей авторегрессии временных рядов при наличии помех», Математическое моделирование, т. 17, вып. 5, сс. 10-16, 2005.
[22] A. E. Permanasari, I. Hidayah, и I. A. Bustoni, «SARIMA (Seasonal ARIMA) implementation on time series to forecast the number of Malaria incidence», в 2013 International Conference on Information Technology and Electrical Engineering (ICITEE), 2013, сс. 203-207.
Technique for Formal Methods of Time Series Forecasting Integration in Data Assimilation
Y.S. Timoshenkova, S.V. Porshnev, N.T. Safiullin
Annotation—The article describes the method developed by the authors for integration of formal methods of time series (TS) forecasting (autoregressive integrated moving average (ARIMA), singular spectrum analysis, group method of data handling, artificial recurrent neural network with long short-term memory) into the Data Assimilation (DA), in cases where mathematical model of the dynamic system generating the TS is not known (for example TS consisting of economic indicators). The performance of the proposed integration method is illustrated on the example of forecasting TS named "Air Passengers" using the DA method based on the ensemble Kalman filter with integrated ARIMA method. Estimations of the accuracy of "Air Passengers" TS forecasts is calculated using ARIMA and the proposed integration method. Article discusses the advantages and disadvantages of the proposed method for integrating formal TS forecasting methods into the DA and also directions for its further improvement.
Key Words—Time series analysis, time series forecast, Data Assimilation.
REFERENCES
[1] G. E. Box, G. M. Jenkins, G. C. Reinsel, i G. M. Ljung, Time series analysis: forecasting and control. John Wiley & Sons, 2015.
[2] L. Nerger i dr., «SANGOMA: Stochastic Assimilation for the Next Generation Ocean Model Applications EU FP7 SPACE-2011-1 project 283580».
[3] M. Kanamitsu, «Description of the NMC global data assimilation and forecast system», Weather and forecasting, t. 4, vyp. 3, ss. 335-342, 1989.
[4] A. V. Kostrov i A. M. Mardjuk, «Ob odnoj matematicheskoj zadache teorii dvizhenija raket», Zhurnal vychislitel'noj matematiki i matematicheskoj fiziki, t. 34, vyp. 10, ss. 1427-1443, 1994.
[5] L. Kilian i M. P. Taylor, «Why is it so difficult to beat the random walk forecast of exchange rates?», Journal of International Economics, t. 60, vyp. 1, ss. 85-107, 2003.
[6] G. Welch i G. Bishop, «An introduction to the Kalman filter», 1995.
[7] G. P. Zhang, «Time series forecasting using a hybrid ARIMA and neural network model», Neurocomputing, t. 50, ss. 159-175, 2003.
[8] N. Golyandina, V. Nekrutkin, i A. A. Zhigljavsky, Analysis of time series structure: SSA and related techniques. CRC press, 2001.
[9] A. G. Ivahnenko, Dolgosrochnoe prognozirovanie i upravlenie slozhnymi sistemami. Tehnika, 1975.
[10] Y. Timoshenkova, S. Porshnev i N. Safiullin, «On the possibility of correction of the forecasting of the Lorenz attractor dynamic characteristics using experimental data and data assimilation», v Journal of Physics: Conference Series, 2018, t. 1053, s. 012004.
[11] Y. Timoshenkova, S. Porshnev, N. Safiullin i O. Ponomareva, «On the Possibility of the Forecast Correction Lorenz Systeme Use Particle Filter», v 2018 International Conference on Applied Mathematics and Computational Science, ICAMCS. NET 2018, 2018, ss. 27-30.
[12] P. L. Houtekamer i H. L. Mitchell, «Data assimilation using an ensemble Kalman filter technique», Monthly Weather Review, t. 126, vyp. 3, ss. 796-811, 1998.
[13] M. Ghil i P. Malanotte-Rizzoli, «Data assimilation in meteorology and oceanography», Advances in geophysics, t. 33, ss. 141-266, 1991.
[14] B. K. W. Lahoz i R. Menard, Data assimilation. Springer, 2010.
[15] E. N. Lorenz, «Deterministic nonperiodic flow», Journal of atmospheric sciences, t. 20, vyp. 2, ss. 130-141, 1963.
[16] M. Tabor, Haos i integriruemost' v nelinejnoj dinamike. 2001.
[17] M. Roth, G. Hendeby, C. Fritsche, i F. Gustafsson, «The Ensemble Kalman filter: a signal processing perspective», EURASIP Journal on Advances in Signal Processing, t. 2017, vyp. 1, ss. 1-16, 2017.
[18] «Air Passengers». https://kaggle.com/abhishekmamidi/air-passengers (prosmotreno dek. 01, 2021).
[19] M. Bocquet, «Introduction to the principles and methods of data assimilation in geosciences», Notes de cours, Ecole des Ponts ParisTech, 2014.
[20] Glavnye komponenty vremennyh rjadov: metod «Gusenica» / Pod red. D.L.Danilova, A.A.Zhigljavskogo. 1997. Prosmotreno: dek. 20, 2021. [Onlajn]. Dostupno na: https://rusneb.ru/catalog/001980_000024_RU_FESSL_MAIN_1331818796 268482813/
[21] A. N. Tyrsin, «Postroenie modelej avtoregressii vremennyh rjadov pri nalichii pomeh», Matematicheskoe modelirovanie, t. 17, vyp. 5, ss. 1016, 2005.
[22] A. E. Permanasari, I. Hidayah, i I. A. Bustoni, «SARIMA (Seasonal ARIMA) implementation on time series to forecast the number of Malaria incidence», v 2013 International Conference on Information Technology and Electrical Engineering (ICITEE), 2013, ss. 203-207.
Yulia Timoshenkova, Ural Federal University, Yekaterinburg, Russia, https://orcid.org/0000-0002-4137-0484 (e-mail: julia.timoshenkova@urfu.ru).
Sergey Porshnev, Ural Federal University, Yekaterinburg, Russia, https://orcid.org/0000-0001-6884-9033 (e-mail: s.v.porshnev@urfu.ru).
Nikolai Safiullin, Ural Federal University, Yekaterinburg, Russia, https://orcid.org/0000-0002-0004-5634 (e-mail: n.t.safiullin@urfu.ru).
The reported study was funded by RFBR according to the research project No 20-37-90006
