CC BY
249
13. Engel, J. Conservatism, accuracy and efficiency: Comparing Value at Risk methods. Discussion paper 2. [Text] /J. Engel, M. Gizicki. -Australia: Australian Prudential Regulation Authority, Reserve Bank of Australia, 1999. - 64 p.
14. Longerstaey, J. RiskMetrics technical document [Text] / J. Longerstaey, M. Spenser. - 4th ed. - New York: J.P. Morgan/Reuters, 1996. - 281 p.
15. Лушавин, А. П. Обработка временных рядов параметров технологических процессов [Текст] / А. П. Лушавин. - LAP LAMBERT Academic Publishing, 2012. - 160 c.
16. Лоскутов, А. Ю. Применение метода локальной аппроксимации для прогноза экономических показателей [Текст] / А. Ю. Лоскутов, Д. И. Журавлев, О. Л. Котляров // Вопросы анализа и управления риском. - 2003. - Т. 1, № 1. - С. 21-31.
17. Чистякова, А. А. Информационная технология прогнозирования нестационарных временных рядов с использованием сингулярного спектрального анализа [Текст] / А. А. Чистякова, Б. В. Шамша // Восточно-Европейский журнал передовых технологий. - 2014. - № 2/4 (68). - С. 24-30. - Режим доступа: \www/URL: http://journals.uran.ua/eejet/article/view/22158.
-□ □-
В роботi запропоновано гiбридт математичн моделi рiзноi складностi для прогнозування часових рядiв i методи iх структурноi iдентифiкацii. Методи прогнозування, що було розроблено, заснован на стль-ному використанн методiв «Гусениця»-ББА та Бокса-Дженктса, а гiбридт математичн моделi, як синте-зовано на iх основi, е прюритетними на сьогодншнш день iмовiрнiсно-детермiнованими декомпозицшними моделями. Експериментальн результати показують ефективтсть запропонованих методiв прогнозування Ключовiслова: прогнозування часовихрядiв, структурна iдентифiкацiя моделi, декомпозицшна модель,
метод Бокса-Дженктса, метод «Гусениця»-ББА □-□
В работе предложены гибридные математические модели различной сложности для прогнозирования временных рядов и методы их структурной идентификации. Разработанные методы прогнозирования основаны на совместном использовании методов «Гусеница»-55А и Бокса-Дженкинса, а синтезированные на их основе гибридные математические модели являются приоритетными на сегодняшний день вероятностно-детерминированными декомпозиционными моделями. Экспериментальные результаты показывают эффективность предложенных методов прогнозирования
Ключевые слова: прогнозирование временных рядов, структурная идентификация модели, декомпозиционная модель, метод Бокса-Дженкинса, метод
«Гусеница»-55А -□ □-
УДК 519.67
|DOI: 10.15587/1729-4061.2014.28172|
гибридные модели и методы прогнозирования временных рядов на основе методов «гусеница^-ssa и бокса-дженкинса
В. Н. Щелкалин
Инженер 1-й категории Кафедра прикладной математики Харьковский национальный университет радиоэлектроники пр. Ленина, 14, г. Харьков, Украина, 61166 E-mail: [email protected]
1. Введение
К настоящему моменту времени создано большое количество математических моделей и методов анализа и прогнозирования временных рядов (ВР) и наметилась тенденция их комбинирования с целью получения лучших характеристик комбинированной модели. А также прослеживается тенденция перевода теоретических наработок одних методов на другие, где это представляется возможным и целесообразным.
В задачах моделирования и прогнозирования ВР различной природы в 70-90-е годы прошлого века получили широкое применение статистические модели. Популярность моделей данного типа объясняется достаточно высокой степенью следующих показателей:
экономность по количеству оцениваемых параметров; высокая скорость процедуры идентификации модели; малая трудоёмкость и ресурсоемкость реализации при использования на доступных тогда ЭВМ малой производительности; наилучшие результаты прогнозирования процессов в различных областях науки и техники в те времена; математическая обоснованность теории; адекватность для решения целого ряда задач теории и практики управления, контроля, моделирования и прогнозирования процессов в различных областях человеческой деятельности. Модель авторегрессии -проинтегрированного скользящего среднего (АРПСС, на англ. - АШМА) является наиболее распространённой из классических статистических моделей прогнозирования.
!© В.
Классические методы прогнозирования обычно предназначены для работы с линейными и стационарными ВР. Однако большинство естественных материальных процессов реальных физических систем являются нелинейными и нестационарными. Поэтому далее наметилась тенденция критического отношения к статистической постановке проблемы идентификации моделей ВР, особенно в случае, когда отсутствует возможность получения представительных выборок для построения математических моделей, статистических характеристик процессов и проверки их адекватности. Кроме того, статистическая теория использует операции осреднения по множеству реализаций, что в целом ряде случаев приводит к ухудшению математической модели, особенно в условиях малых и нестационарных выборок [1, 2]. В этих случаях эффективно использовать алгебраический, детерминированный, а не статистический подход к решению проблемы идентификации. И лишь в последние десятилетия достаточно активно начали развиваться методы прогнозирования линейных, но нестационарных ВР и нелинейных, но стационарных ВР. И относительно недавно начали появляться гибридные модели прогнозирования, синтезируемые достаточно сложными комбинациями методов, учитывающих как детерминированные, так и стохастические составляющие процессов.
Таким образом, существуют три основные группы моделей и методов моделирования ВР:
1) детерминированные (алгебраические);
2) вероятностные (статистические);
3) комбинированные вероятностно-детерминированные.
Актуальным представляется проведение синтеза классов математических моделей, основанных на совместном использовании идей двух достаточно хорошо теоретически обоснованных методов: детерминированного метода «Гусеница»-$$А и статистического метода Бокса-Дженкинса. Это связано с очевидными преимуществами метода «Гусеница»-SSA, среди которых можно выделить следующие: данным методом можно разложить ВР на интерпретируемые аддитивные составляющие, не требуя при этом стационарности ряда и знания модели тренда, а также сведений о наличии в ряде периодических составляющих и их периодах. Модели метода Бокса-Дженкинса включают в себе все приведенные выше преимущества статистических моделей прогнозирования. Однако данные методы имеют и ряд недостатков. Так, например, недостатками метода Бокса-Дженкинса можно считать заниженное количество оцениваемых параметров, используемых в модели, искажение долгосрочной зависимости ВР и недостаточное количество публикаций со сравнительным анализом прогнозирования данным методом. Среди недостатков же метода «Гусеница»-SSA можно выделить использование неоптимального с точки зрения точности воспроизведения некоторых ВР ортогонального базиса векторов траекторной матрицы ВР, существенно неавтоматическая группировка компонент сингулярного разложения траекторной матрицы ряда для получения составляющих исходного ряда, отсутствие модели не позволяет проверять гипотезы о наличии в ряде той или иной составляющей. Данный метод является непараметрическим. Для проверки подобных гипотез требуется построение модели, ко-
торое в свою очередь, может быть проведено на основе информации, получаемой с помощью метода «Гусени-ца»-SSA. Этот непараметрический метод позволяет получить результаты, часто лишь незначительно менее точные, чем многие параметрические методы при анализе ряда с известной моделью. Совместное же использование данных методов приводит к синергии, взаимно компенсируя противоположные по природе недостатки входящих в него моделей.
В статье приведены модели трендового подхода к прогнозированию ВР. Трендовый подход заключается в моделировании процесса как отклонения фактических значений ВР относительно трендовой составляющей, в роли которой в предлагаемой модели выступает линейная рекуррентная формула (ЛРФ) метода «Гусе-ница»-SSA.
Однако реальные ВР в различных отраслях науки и техники являются сложными и обычно порождаются двумя или более процессами, которые могут изменяться во времени. Каждый из этих процессов может быть описан различными моделями. Качественными признаками сложности являются нелинейность процесса, неопределённость его математической модели, наличие внешних возмущений. В результате перечисленных свойств построенная традиционными методами математическая модель сложного процесса, как единого и неделимого, является практически не реализуемой и малопригодной для использования в задачах прогноза. Поэтому также в данной статье предлагаются модели и методы декомпозиционного подхода к прогнозированию, основанные на использовании методов «Гусеница»-SSA и Бокса-Дженкинса. Основной целью декомпозиционного подхода к прогнозированю является разделение исходного ВР на множество ВР, с более простой структурой, рассматриваемых независимо друг от друга. При этом декомпозиционная модель строится следующим образом:
- осуществление декомпозиции, т. е. непосредственно разложение исходного ВР;
- построение отдельных упрощённых, как правило, независимых моделей, соответствующих компонентам разложения;
- реализация общей декомпозиционной модели процесса как объединения построенных упрощённых моделей.
В статье производится сравнительный анализ предложенных математических моделей и методов прогнозирования на ВР потребления электроэнергии и природного газа.
2. Анализ литературных данных и постановка проблемы
В рассмотренной литературе [3-34] многочисленные исследования направлены на анализ математических моделей и методов среднесрочного и долгосрочного прогнозирования ВР в различных областях науки и техники. В [2] представлены аддитивные, мультипликативные и смешанные модели прогнозирования электропотребления. В [1] рассмотрены различные комбинированные математические модели для прогнозирования трендовой, сезонной, недельной, остаточной составляющих процессов потребления элек-
троэнергии и предложены декомпозиционным метод моделирования (ДММ) и метод окаймляющих функций. Среди предлагаемых моделей прогнозирования отдельных составляющих процесса используются экспоненциальное сглаживание, модель Хольта-Винтер-са, двойного экспоненциального сглаживания, фильтр Калмана, ряды Фурье, модели спектрального анализа, экспоненциально взвешенного скользящего среднего и др. Для оценивания эффективности полученных прогнозов модели обычно сравнивались с моделями ARIMA [4-7].
ВР потребления целевых продуктов в системах энергетики являются существенно нелинейными со средне- и долгосрочными зависимостями, поэтому аппроксимация таких сложных ВР только моделью ARIMA не всегда является удовлетворительной. Было доказано, что для повышения точности прогнозирования нестационарных ВР необходимо эффективно извлечь трендовую, гармонические и шумовую составляющие и далее подставлять их в модели прогнозирования. Таким образом, методы декомпозиции ВР являются благоприятными для первичной обработки данных.
Сингулярный спектральный анализ (SSA) является одним из наиболее эффективных методов для анализа ВР со сложными периодическими компонентами. Зарождение SSA обычно ассоциируется с публикацией работы [8], а также метод SSA независимо разработан в России и в нескольких группах в Великобритании и США. SSA стал широко используемым методом анализа периодических и трендовых составляющих ВР [9-12]. В последнее время были разработаны новые модели для прогнозирования, основанные на SSA. В [13] разработана модель на основе сингулярного спектрального анализа и стандартной искусственной нейронной сети. В [14] создана гибридная модель прогнозирования ВР на основе метода опорных векторов (SVM) в сочетании с SSA. В [15] предложили линейную рекуррентную формулу сингулярного спектрального анализа (SSA-LRF), и эта модель была использована для решения некоторых практических задач [16, 17].
В последнее время появилось достаточно большое количество работ, посвященых и другим гибридным декомпозиционным моделям и методам прогнозирования. В [18] SVM комбинировалась с вейвлет-анализом. В [19] предложен двухэтапный адаптивный подход для прогнозирования ВР потребления электроэнергии. В этом подходе на первом этапе методом EMD (Empirical Mode Decomposition) производится декомпозиция прогнозируемого ВР на собственные модальные функции и применяется к ним преобразование Гильберта. На втором этапе, полученные модальные функции и их мгновенные амплитуды подаются на вход прогнозной нейронной сети. В [20] разработан комбинированный метод прогнозирования тоже основанный на разложении Хуанга и искусственной нейронной сети. В [21] был предложен подобный метод, но в качестве метода разложения использовалось разложение на независимые компоненты ICA (independent component analysis). В [22] комбинируются модели SVR (support vector regression) и ARIMA. Достаточно много статей посвящено гибридным моделям, сочетающим нейронные сети с моделями ARIMA [23-26]. Гибридные модели
на основе SVM и ARIMA рассмотрены в [27-29]. В [30] EMD комбинируется с моделями ARIMA. В [31] ARIMA комбинируется с вейвлетами.
В данной статье разрабатываются гибридные математические модели и методы на основе методов «Гусеница»-SSA и Бокса-Дженкинса. В [32] предложена комбинированная аддитивная модель прогнозирования, в которой к ЛРФ добавлялась модель авторегрессии - проинтегрированного скользящего среднего. Исследования комбинирующие метод «Гусеница»-SSA с другими моделями можно встретить и в зарубежной литературе. Так, например, в [33] была разработана подобная модель, в которой ЛРФ комбинировалась с моделью авторегрессии и применялась такая гибридная модель для краткосрочного прогнозирования нагрузки иранской национальной энергосистемы. Другой примером совместного использования метода «Гусеница»-SSA с моделями Бокса-Дженкинса можно найти в [3]. Там методом «Гусеница-SSA производилось разложение исходного ВР, затем производилась группировка ВР разложения на интерпретируемые составляющие, такие как трендовая, сезонные и шумовые. Далее для каждого сгруппированного ВР идентифицировалась и оценивалась своя модель ARIMA. Остаточная ошибка моделирования корректировалась авторегрессионной моделью и находились прогнозы синтезированной моделью. Примерно в то же время в [34] была предложена подобная модель, однако, в отличие от [3], для идентификации и оценивания сгруппированных ВР использовались сезонные модели авторегрессии - проинтегрированного скользящего среднего (САРПСС, на англ. - SARIMA), а также такая гибридная модель была обобщена уже для случая прогнозирования ВР, зависящего от нескольких экзогенных ВР.
Таким образом, анализ рассмотренной литературы позволяет сделать вывод, что в настоящее время происходит отход от статистической постановки задачи прогнозирования и эффективное прогнозирование нестационарных ВР предполагает использование различных приёмов декомпозиции, синтеза прогнозных моделей и отбора их из множества альтернатив.
3. Цели и задачи исследования
Целью проведенных исследований является разработка перспективных по мнению автора декомпозиционных вероятностно-детерминированных математических моделей и методов их идентификации для прогнозирования нестационарных ВР со сложной структурой и проверка эффективности прогнозирования разработанными моделями в сравнении с прогнозами, получаемыми вероятностными моделями SARIMA, обобщёнными для прогнозирования ВР с несколькими сезонными составляющими [4].
Для этих целей в работе использованы методы «Гусеница»-SSA и Бокса-Дженкинса и на их основе разработаны гибридные математические модели и методы прогнозирования нестационарных ВР. Метод «Гусеница»-SSA применялся для разложения исходного ВР на ВР с более простой структурой, а метод Бокса-Дженкинса - для построения математической
модели данных компонент разложения. Выбор различных параметров метода «Гусеница»-$$А приводил к созданию в сочетании с методом Бокса-Дженкинса различных вариантов методов прогнозирования и к синтезу на их основе гибридных математических моделей с различными структурами.
4. Математические модели и методы прогнозирования нестационарных временных рядов
В данном разделе сначала приведём короткое изложение двух классических методов прогнозирования нестационарных ВР: детерминированного метода «Гу-сеница»-SSA и статистического метода Бокса-Джен-кинса. Приведенные в них сезонная модель авторегрессии - проинтегрированного скользящего среднего (САРПСС) метода Бокса-Дженкинса и линейная рекуррентная формула (ЛРФ) метода «Гусеница»-SSA необходимы далее для описания гибридных моделей на их основе.
4. 1. Классические методы прогнозирования нестационарных временных рядов
4. 1. 1. Метод «Гусеница»-88А
Рассмотрим ВР yt =(у1,у1,—,уп) длины п. Основным параметром метода «Гусеница»-SSA является параметр L, 1<Ь<п и называется длиной окна. Базовый алгоритм метода «Гусеница»-SSA состоит из следующих этапов [35].
1. Разложение ВР.
1. 1. Вложение. На начальном этапе исходный ВР у1, 1 = 1,11 переводится в последовательность К=п-Ь+1 многомерных векторов вложения размерности Ь:
Ч =(у„у,-1 ...,у1+^1 )т, 1-1-К.
Затем из них строится так называемая траекторная матрица Y следующим образом:
Y = [У[ Y2 ... Yк ]
Y=(уй
^=1
УО У1 У 2 . •• УК-1
У1 у 2 у3 . •• УК
У 2 у3 у 4 . •• ук+1
у L-1 у L у L+1
уп-
1
1. 2. Сингулярное разложение. К траекторной матрице ВР применяется SVD-разложение:
Y = = Е Y1, (1)
1=1 1=1
где Yi = ^"и^Т ; и1, и V , 1 = - левые и правые сингулярные вектора, соответствующие упорядоченным по убыванию собственным числам ^ , 1 = матрицы 8 = YYT , d=rank Y .
2. Восстановление ВР.
2. 1. Группировка. Всё множество индексов {1, 2, ..., Ц делится на г непересекающихся подмножеств 11, 12, ..., 1Г. Для каждого набора индексов 1к = {1(к)д2к),...д[,к1)} ,
к = 1,г определяется матрица YI = ЕУ(к) . Тогда разло-
к ]=1 1
жение (1) можно записать следующим образом:
У = Е У1к.
Процедура выбора множеств 11, 12,...,1Г и называется группировкой собственных троек.
2. 2. Диагональное усреднение. На этом этапе каждая матрица У1, к = 1,г переводится во ВР g[k) = (g(k),g2k),.,gnk)), сопоставляя каждой к-ой «диагонали» матрицы среднее арифметическое значений данной «диагонали».
Применив диагональное усреднение к матрицам У1ь , к = 1,г получаем ряды у[к' = (у1к),у2к),.--,уПк)), к = 1,г. Тогда исходный ряд раскладывается в сумму г рядов:
Уt
-у t = 1 у 1°.
(2)
3. Прогноз. _
Обозначим wi = у ¡, 1 = 1,п. Рекуррентный прогноз на 1 шагов вперед определяется по следующей формуле:
L-1 _
w1 = Е gjWi_j,i = п + 1,п +1, (3)
]=1
где коэффициенты gj, j = 1^-1 являются элементами вектора = ) , определяемого по формуле:
1 -V
■е^иг,
где п - последняя компонента вектора и1, 1 = 1,г; иг, 1 = 1,г - вектор, состоящий из первых L - 1 компонент вектора и1; V2 = п2 + тс2 + ... + ТС;2. Формула (3) называется линейной рекуррентной формулой (ЛРФ).
Имеется и векторный вариант прогнозирования методом «Гусеница»-SSA, который в ряде случае может давать более точные прогнозы в сравнении с вариантом метода рекуррентного прогнозирования [35].
4. 1. 2. Модель авторегрессии - проинтегрированного скользящего среднего
В методе Бокса-Дженкинса уделено особое внимание проблеме выбора модели и её оцениванию. В этом методе используются идеи, что нестационарные ВР можно преобразовать в стационарные путём перехода от исходного ВР к его разностям соответствующего порядка d=1, 2, .... Преобразованный таким образом ВР можно далее описать моделями авторегрессии (АР), скользящего среднего (СС) или смешанной АРСС моделью соответствующих порядков.
Сезонная модель АРПСС (САРПСС, на англ. -SARIMA) в операторной форме может быть представлена в следующем виде [4]:
Уt =
е;.(В) ф++(ВГ
(4)
где у^ 1 = 1,п исходный или преобразованный (нормированный или прологарифмированный) центрированный ВР; п - объём выборки; В - опе-
1
ратор сдвига по времени на одну единицу назад, такой что В'хк = хк-1; ( ) - обобщенный оператор
п. nS
авторегрессии порядка ;+= + ^DiSi, ;* = ^piSi;
1=1 1=1
Ф++(В) = Ф;*(В)У5С/У5С22 D1, 1 = - порядок
взятия разности Si; Si, 1 = 1,п5 - период 1-й периодической компоненты, причем S1 = 1; п5 - количество периодических компонент; У5 и В5 - упрощающие операторы такие, что У5 х1 =(1 -В51)х1 = х1 -х1-5.
Ф*. (В) - обобщенный оператор авторегрессии порядка ;* вида Ф*;.(В) = ПФ';1 (В51); Фр (в5), 1 = Т^П^ -
полиномы от В5' степеней ;1 соответственно, определяющие составляющие авторегрессии периодических компонент с периодами соответственно;
0*. (В) - обобщенный оператор скользящего сред-
П5 П5 . .
него порядка = ^вида 0*. (В) = ^0* (В5);
1=1 * 1=0 1
0* (В5 ), 1 = 1,п5 - полиномы от В5 степеней *1 соответственно, определяющие составляющие скользящего среднего периодических компонент с периодами 51 соответственно;
а1 - случайный процесс типа белый шум. Идентификация параметров ;1, 1 = 1,п5 осуществляется путём анализа частной автокорреляционной функции (ЧАКФ), а параметров *1, 1 = 1,п5 - путём анализа автокорреляционной функции (АКФ). Параметры п5, Di, также идентифицируются путём анализа АКФ и ЧАКФ процесса. Автоматизировать процесс структурной идентификации таких моделей можно перебором моделей с различной структурой и отбором наилучших по информационному критерию Акаике ^Ю) и байесовскому информационному критерию Шварца (В1С).
4. 1. 2. 1. Алгоритм обобщённого метода Бокса-Дженкинса
Блок-схема алгоритма, обобщённого на случай нескольких сезонных составляющих ВР, метода Бок-са-Дженкинса (рис. 1) состоит из следующих этапов [4]:
1) начало;
2) ввод обучающей выборки значений исходного ВР yí в дискретные моменты времени 1 = 1,п;
3) идентификация порядков взятия разностей и преобразование ряда. На основании анализа вида автокорреляционной функции процесса yí определяем порядок взятия разности Dn с периодом 5 =1. Определяем периоды периодических компонент 51 с убывающими амплитудами. На основании анализа автокорреляционных функций процессов
у'=У5с;у1; у'Ь=У5С; у£уГ. у1!=у5С/ у5С22-уЗ^,
определяем порядки взятия разностей , 1 = 1,п5 -1. По известным значениям , преобразуем исходную выборку yí в
^=у5с; у5с; х- х^,
4) идентификация параметров модели ;1, *1, 1 = 1,п5 . Осуществляется по виду автокорреляционной и частной автокорреляционной или автоматически перебором моделей с различной структурой и отбором наилучших по информационному критерию Акаике ^Ю) и байесовскому критерию Шварца (В1С);
5) грубая оценка параметров модели. Как показали экспериментальные исследования сходимости алгоритма параметрической оптимизации, для правильно идентифицированной модели сходимость этого алгоритма обеспечивается при задании в качестве грубых оценок параметров модели констант типа ±0.1;
6) точная оценка параметров модели. Осуществляется в соответствии с алгоритмом параметрической оптимизации. В данной работе использовался метод Левенберга-Марквардта. В качестве начального приближения итерационной процедуры вначале используются грубые оценки параметров модели, в дальнейшем - точные оценки параметров модели, полученные для предыдущей обучающей выборки;
7) проверка адекватности полученной модели. Осуществляется в соответствии с критерием
а М£ Г„ (к) ,
(5)
где гаа (к) - корреляционная функция остаточных ошибок а1 модели.
Сравнив значение данного критерия с табличными значениями %2 при заданном уровне значимости, приближенно проверяется гипотеза об адекватности полученной модели. Если полученная модель неадекватна исходным данным, то осуществляется переход к четвертому этапу. Если модель адекватна - к следующему этапу;
8) вычисление прогноза у1 (1), 1 = с заданным упреждением L. Осуществляется по разностному уравнению
у<(1) = 0о + 1Ф+ Еу 1+1-1 ]-£0*[а1+ч] ,
1=1 ]=1
где 00 - общая константа, позволяющая учесть детерминированный тренд. Условные математические ожидания величин у и а в моменты време-в ни 1+1-1 и +1- соответственно равны
[уы-1 ] =
[а1+Ч] =
у! (1 -1),1 > 1;
у1+1-1,1 * 1,
0,1 >
у1+ч- Уы-ц(1),1 * *
где у1 (1 -1) - прогноз величины у1, вычисленной в момент времени 1 с упреждением (1 -1); Ф+, 1 = 1,; + и 0,, ] = 1,** - коэффициенты полиномов Ф;+ (В) и 0*. (В) соответственно.
Дисперсия прогноза V (1), 1 = , а также его вероятностные пределы находим соответственно по уравнениям
V (1)=о2 е V,2,
)=0
у,+1 (±) = (1)± и^,
где ие/2 - квантиль уровня 1 -е/2 стандартного нормального распределения. Веса V; определяются по формуле
упреждением в единицу. Здесь же выполняем необходимые преобразования новой обучающей выборки уы в
I 3, у Б,
V. = 1
1,] = 1;
е -а;,] > 1
¡=1
1>Р+; й* = о,
где Ф + = 0 ,
9) печать прогноза у,.(1), 1 = 1,Ь; "
10) ввод фактического значения уы , ¿ = 1,Ь ;
И) вычисление остаточной ошибки. При вводе в систему каждого нового значения уы, 1 = 1, Г- определяется остаточная ошибка модели аы = уы У (1)
12) проверка условия: если /<1, то осуществляется переход к следующему этапу, если нет - к 15-му этапу;
13) коррекция прогнозов. Используя веса и остаточную ошибку модели аы в момент времени 1+1, производим коррекцию прогнозов уДк), к = 1 + 1,Ь в соответствии с уравнением
УмОЬ^М + т+Г.
14) печать скорректированных прогнозов и возврат к 10-му этапу;
15) проверка условия продолжения вычислений. Если условие выполняется, то осуществляется переход к следующему этапу, если нет - к 17-му;
16) формирование новой обучающей выборки и нового ряда остаточных ошибок модели. Осуществляется путем отбрасывания Ь первых и прибавления Ь новых реализаций процесса и фактических ошибок прогнозов с
17) проверка адекватности модели.
Для ряда остаточных ошибок модели, соответствующих новой обучающей выборке, вычисляем значение критерия (5) и проверяем адекватность модели. Если модель адекватна, переходим к 8-му этапу. В противном случае, предполагая, что структура модели не изменилась, а изменились только значения параметров модели, переходим к 6-му этапу. Если новые данные привели к изменению стохастической структуры модели, это будет обнаружено на 7-м этапе (рис. 1).
Рис. 1. Блок-схема алгоритма обобщённого метода Бокса-Дженкинса
Этот алгоритм обеспечил высокую точность прогнозов процессов потребления целевых продуктов в системах энергетики такими потребителями, для которых превалирующее влияние оказывают хронологические факторы [4].
4. 2. Гибридные математические модели и методы прогнозирования нестационарных временных рядов на основе методов «Гусеница»-88А и Бокса-Дженкинса
В данном разделе рассмотрены трендовый и декомпозиционный подходы к прогнозированию ВР. Для каждого из подходов предложены соответствующие математические модели и методы их построения, основанные на совместном использовании двух методов: детерминированного метода «Гусеница»-SSA и статистического метода Бокса-Дженкинса.
4. 2. 1. Трендовый подход к прогнозированию нестационарных временных рядов
В общем случае ВР у1 (1=1, 2,..., п) может быть разложен на трендовую составляющую уТ (включая её разновидности - сезонную или гармоническую) и случайную составляющую е^
у1 = уТ+е1.
Очевидно, необходимо прогнозировать как трендо-вую, так и случайную составляющие.
Трендовый подход состоит в моделировании ВР как отклонения фактических значений ВР относительно его трендовой составляющей. Для прогнозирования трендовой составляющей наиболее широкое распространение получили следующие методы: метод наименьших квадратов (МНК), взвешенный МНК, аппроксимация полиномами, метод гармонических весов, метод экспоненциального сглаживания, адаптивного сглаживания, метод построения сплайн-функций, компенсация трендовой составляющей на основе вычислений разностей п-го порядка, методы скользящего среднего, методы оптимальной фильтрации, методы моделирования конечными рядами Фурье, спектральные методы, методы декомпозиции по ортогональным векторам или функциям, методы имитационного моделирования, методы нейро-сетевого и нечёткого моделирования, методы вариационного исчисления, метод группового учёта аргументов, метод «Гусеница»-SSA [1, 2] и др. Для прогнозирования же случайной составляющей используются методы вероятностного моделирования, регрессионные статистические модели, метод Цепей Маркова, модели экспоненциального сглаживания Брауна, взвешенного среднего, Хольта-Уинтерса, обобщённый фильтр Винера (метод прогнозирования Заде-Рагаззини), спектральные статистические модели (канонические разложения, разложение Карунена-Лоэ-ва, модели компонентного и факторного анализов), прогнозирующий фильтр Калмана-Бьюси, модели кластерного анализа и теории распознавания образов [1, 2] и др.
4. 2. 1. 1. Гибридные математические модели трен-дового подхода к прогнозированию нестационарных временных рядов на основе методов «Гусеница»-88А и Бокса-Дженкинса
В настоящее время всё большее распространение начинают получать математические модели прогнози-
рования, являющиеся комбинацией статистических и детерминированных моделей.
Вернёмся к подпункту 1 этапа восстановления ВР метода «Гусеница»-SSA возьмём параметр г=1, т. е. из множества индексов {1, 2,., L} отберём некоторые или все индексы в подмножество 1Т с {1, 2,., L}. Тогда на этапе диагонального усреднения метода «Гусени-ца»-SSA исходному ВР будет поставлен в соответствие некоторый отфильтрованный (сглаженный) ряд уТ. Т. е., согласно (2), имеем у1 = уТ . Вычислим остаточный ВР у а как у а = у1 - у Таким образом, исходный ВР будет разложен на следующую сумму:
у1=у Т+у а.
Данное разложение нужно стараться производить так, чтобы ВР у^ соответствовал тренд-сезонной составляющей исходного ВР, а остаточный ВР у а - шумовой. Это производится путём анализа собственных значений и собственных векторов траекторной матрицы [35]. При этом могут возникнуть трудности со слабой разделимостью некоторых ВР [35]. _
Вспомним про обозначение = уТ, 1 = 1,п и, для удобства записи, введём обозначение: а1 = у а.
В основе предлагаемой математической модели трендового подхода к прогнозированию ВР, основанного на совместном использовании методов «Гусени-ца»-SSA и Бокса-Дженкинса, лежит идея добавления к ЛРФ, вычисляемой методом «Гусеница»-SSA, модели АРПСС, определяемой методом Бокса-Дженкинса, построенной для остаточных ошибках а1 модели, определяемых путём вычитания из исходного ВР значе-
ний Введём обозначение g(В) = ^giBi, где
полином от введенного выше оператора задержки В, коэффициенты gi, 1 = -1 которого равны коэффициентам ЛРФ = ^giwt-i метода «Гусеница»-SSA.
1=1
Если тренд-сезонную компоненту и остаточные ошибки прогнозировать отдельно, тогда описанная выше модель будет иметь вид (6), иначе - вид (7):
е1 = у1 - wt,t=1,п;
е, =
0* (в) ,
УФ; (ВУ
У1 (1) = - (1К (1),
(6)
где ф;(В) = (1 -ф1В-ф2В2- — -ф;В;) - оператор авторегрессии порядка ;; 0Ч (В) = (1-01В -02В2 -- -0ЧВЧ) -оператор скользящего среднего порядка *; у1 (1), (1) и е1 (1) - прогнозы ВР у1, и е1 соответственно в момент времени 1+1. Прогнозы (1) и (1) вычисляются путём перехода от операторных уравнений
с 0* (В) 1 УФ; (В)"1
а и — =
к разностным:
et (l) = e0 + (ф4 + l)[et+,-1 ] +
+ (ф2-ф1 )[£t+i-2 ] + (ф3 ф2 )Х
x[et+i-3 ] + ••• + (фр -фр-1 )[еы_р ] +
+фр [et+i-(p+D]-^le. [at+i-i ],
Yt =
Yt =
+ e;-(B) ^ +— / ч a, ,
t ФР+(В) t,
e;.(B) ' ф'++(в)at,
Yt"
(10)
(11)
где
Yt (l - i),l > i
Yt+i-i, l * i, 0,l > j,
Yt+i-j - Yt+i-j-1 (1),l *j
wt (i)=z gi [wt+i-i],
где
wt (i - i),i >i,
[wt+i-i ]=l
Lwt+i-i,i *i.
Запишем более общую модель:
Yt =
MB! a
Уфр (B) t
(7)
где параметр L - длина окна метода «Гусеница»-SSA; ф; (В) - оператор авторегрессии порядка ;; 0* (В) -оператор скользящего среднего порядка *; а1 - случайный процесс типа белый шум.
Оператор g (В) может применяться непосредственно к процессу у1 , тогда коэффициенты полиномов 0* (В) и ф; (В) и ВР а1 изменятся, а модель запишется следующим образом:
Yt = g(B) Yt"
eq (B) a'
уФ',(в) t
(8)
Если теперь ввести обозначение g,(B) = 1 -модель (8) примет вид:
g'(B) Yt =
eq (B) ,
уф'. (BY
(9)
Поясним наличие в формулах (7)-(9) оператора У. После применения метода «Гусеница»-SSA остаточный ВР либо сразу становится стационарным, либо - после однократного применения оператора У . Многократное же применение оператора У может привести к появлению кратных корней характеристического полинома g,(B) =0 модели, лежащих на границе устойчивости. Метод «Гусеница»-SSA в такой модели может применяться для выявления скрытых периодичностей во ВР.
Обобщим модели (7)-(9) добавляя к ЛРФ модель САРПСС, определяемую обобщённым методом Бок-са-Дженкинса как в разделе 4. 1. 2. 1:
g'(B) Yt =
e;-(B) ф 'p+(B)a
(12)
Назовём модель (10) моделью ЛРФ-SSA + SARIMA. Метод «Гусеница^-SSA описывает достаточно широкий класс ВР, представляющих собой сумму произведений полиномов, экспонент и гармоник, подклассом
которых являются ВР вида wn = ^Pke"k" cos(2nrnkn + фк).
к
Поэтому выделенная из исследуемого процесса методом «Гусеница^-SSA детерминированная составляющая wt исходного ВР достаточно хорошо описывается моделями САРПСС. Следовательно, целесообразно рассмотреть следующую модель:
eWw (B) w
wt =—3—a, ;
t ^ww (B) t
et = Yt- wt;
£t = v<!p(i3)V; Yt (i) = w (i) + £t (i).
или модель
=jw(B) e4 (B)
^ (B) at ' УФр (B)
at,
где а- , а1 - случайные процессы типа белый шум (остаточные ошибки моделей). Или в обобщённом виде эти модели запишутся следующим образом:
0С (В)
1 (В) 1 ,
е1 = у1 - wt, „_0**(В) а
t_ фр+н
Yt (i) = wvt (i)+et (i)
ew-(B) w, e;-(B) ф ww++(B)at ф++Нat.
Yt =
(13)
(14)
Дадим модели (14) обозначение Тренд-SSA-SARIMA и поясним. В этой модели, в отличие от (10), тренд-сезонная составляющая описывается не ЛРФ, а моделью SARIMA.
at+i-j =
и
и
4. 2. 1. 2 Метод на основе совместного использования методов «Гусеница^-SSA и Бокса-Дженкинса трендового подхода к прогнозированию нестационарных временных рядов
Ниже приведен алгоритм метода на основе совместного использования методов «Гусеница^-SSA и Бокса-Дженкинса трендового подхода к прогнозированию нестационарных ВР, а также его схема (рис. 2).
1. Применяя этапы 1 и 2 метода «Гусеница^-SSA к исходному ВР, определяется ВР wt = yt , интерпретируемый как тренд-сезонная составляющая исходного ВР.
2. Для wt определяется ЛРФ, управляющая данным рядом, если используем модели (10)-(12) или производится структурная и параметрическая идентификация модели САРПСС (используя алгоритм раз-
at , если же исполь-
е"; (В)
дела 4. 1. 2. 1.) для ВР =фЦВ)<
зуется модель (13) или (14).
3. Вычисляется остаточный ВР е1 следующим образом:
е 1 = у1 - , если используется модель (10);
е 1 = у1 - g(В)Уt для моделей (11), (12);
е 1 = у1 - , если же используется модель (13) или
(14).
4. Строится САРПСС модель остаточного ВР е1:
е\ (В)
е1 =—*+ . . ае (используя алгоритм раздела 4. 1. 2. 1). Ф++(В)
5. Для и е1 вычисляются прогнозы (1) и е1 (1) с заданным упреждением 1 для модели (13), используя соответствующие модели САРПСС.
6. Синтезируется общий прогноз у 1 (1) = (1) + е1 (1) как сумма частных прогнозов для модели (13). Или вычисляется прогноз по модели (14), предварительно перейдя от операторного к разностному уравнению.
4. 2. 2 Декомпозиционный подход к прогнозированию нестационарных временных рядов
Под термином «декомпозиция» при моделировании сложных ВР понимают разложение исходного ВР у1 на более простые у1, 1 = 1^, но при этом совокупность этих более простых процессов эквивалентна исходному процессу:
N
у1 = е у1.
1=1
Рассмотрим декомпозиционный подход в контексте концепции стандартизованного моделирования
ВР, которая заключается в моделировании ВР у1 как
т
совокупности детерминированного тренда у1 , различных колебательных составляющих у3 , 1 = 1,п3 , модулированных по амплитуде и частоте и остаточной составляющей е:
у1 = ут + е у?' + е1.
1=1
Декомпозиционный подход, как правило, используют при моделировании сложных ВР, когда использование традиционных подходов резко затруднено. В настоящее время эффективное моделирование процессов подразумевает использование различных приёмов декомпозиции модели, что позволяет повысить точность и адекватность моделирования в случае нелинейных нестационарных процессов, упростить и повысить устойчивость процесса идентификации [1].
На практике наблюдаемые ВР являются агрегированным представлением нескольких процессов, что ограничивает эффективность использования классических статистических моделей и методов для их анализа и прогнозирования. В большинстве случаев ВР, зависимые от этих порождающие процессов, могут быть разделены с применением гребенки фильтров, методов спектрального анализа, восстанавливающих определенную детерминированную часть ВР.
Исходный временной ряд
П г-
Выделение методом SSA тренд-сезонной составляющей исходного в. р.
1/ш\/ w' V--1 \ J
Вычисление остаточной состовляющей исходного в. р.
Определение ЛРФ для моделей (10)—(12) или SARIMA для моделей (13)—(14)
Вычисление модели »SARIMA для остаточной составляющей в. р.
Синтез модели ЛРФ-SSA+SARIMA
Проверка адекватности синтезированной модели
Вычисление прогноза синтезированной моделью
Рис. 2. Схема трендового подхода к прогнозированию нестационарных временных рядов на основе совместного использования методов «Гусеница^-SSA и Бокса-Дженкинса
Затем отфильтрованный ВР может быть спрогнозирован моделью, агрегированной при помощи специальных процедур отбора.
Таким образом, сложность процессов ассоциируется со следующим:
- значительное влияние внешних факторов;
- нелинейность, нестационарность или априорная неопределённость динамического поведения процесса.
4. 2. 2. 1 Гибридные математические модели декомпозиционного подхода к прогнозированию временных рядов на основе методов «Гусеница»-88А и Бокса-Дженкинса
Для построения адекватных моделей нелинейных и нестационарных ВР необходимо иметь возможность формировать адаптивный базис, функционально зависящий от содержания ВР. Такой подход и реализуется в предлагаемом декомпозиционном подходе к прогнозированию на основе методов «Гусеница»-SSA и Бокса-Дженкинса.
Метод «Гусеница»-SSA используется для поиска скрытых взаимосвязей внутри ВР, являющегося смесью нескольких ВР, не имея никаких априорных знаний
0 механизме их смешивания. Для этого исходный ВР разлагается на несколько ВР с более простой структурой, применением к нему метода «Гусеница»-SSA. Далее каждый из ВР разложения прогнозируется методом Бок-са-Дженкинса. Затем вычисляются прогнозы исходного ВР путем агрегирования прогнозов ВР разложения.
В подпункте 1 этапа восстановления ВР метода «Гу-сеница»-SSA возьмём параметр г=Ь, т. е. всё множество индексов {1, 2,., L} разделим на L подмножеств по одному элементу в каждом следующим образом 11 = {1} ,
1 = . Тогда на этапе диагонального усреднения метода «Гусеница»-SSA исходный ВР будет разложен на следующую сумму:
у1 = ХУ 11).
1=1
Далее декомпозиционный подход к прогнозированию нестационарных ВР на основе совместного использования методов «Гусеница»-SSA и Бокса-Джен-кинса предполагает идентифицировать и строить модели САРПСС каждого ряда разложения у^ , 1 = 1,Ь.
Таким образом, гибридная модель декомпозиционного подхода к прогнозированию ВР на основе методов «Гусеница»-SSA и Бокса-Дженкинса можно записать в следующем виде:
Теперь рассмотрим данную модель в контексте концепции стандартизованного моделирования нестационарных ВР. В подпункте 1 этап восстановления ВР метода «Гусеница»-SSA оставим без изменения как это описано в разделе 3. 1. 1. Т. е. Всё множество индексов {1, 2,., L} делится на г непересекающихся подмножеств 11, 12,..., 1г. Тогда на этапе диагонального усреднения метода «Гусеница»-SSA исходный ВР будет разложен на следующую сумму:
у1 = 1 у 1° ,
(17)
если в непересекающиеся подмножества 11, 12,...,1г вошли все элементы из множества {1, 2,., L}. Или
у1 у 1
(1)
(18)
в противном случае. Причём ВР у'1), 1 = 1,г могут интерпретироваться как трендовые, квазипериодические или шумовые составляющие исходного ВР.
Тогда гибридная модель прогнозирования нестационарных ВР, основанная на концепции стандартизованного моделирования и на методах «Гусеница»-SSA и Бокса-Дженкинса для случая (17) примет вид
у 1°=-
0у
Ф
*у(1), (В) й -
*-ау ,1 = 1,г ,
.А
у1=х у 11),
(В)
У1 (1)=х у 1°(1).
Или
у1 = хтуз
' 0уЙ*у(1)* (В)
1 Фу +уй+(В)
1у 1
(19)
(20)
А для случая (18) -
у ауй,1 = й
у 1 Фу +уй+(В) 1
у 1 = х у 1°;
0у
.(в) , -
(1)= ^ 7 ауй,1 = 1,L ,
У1 =
Фу
(В)
у1 = х у 11),
е1 = у1 - у t,
0**(В) Ф++(В)
у! (1)=Х У 11)(1).
Или
(15) У1 (1) = Х У^К (1).
1=1
Или
(21)
Ь 0уЙ*уй*(В) М
ау .
у( = х фу(1>РуВ+ (В)
(16)
у( = х
0уЙ*уй*(В) й 0** (В)
*у_ау11) + * ' а
х ФуЙ+уй+(В) 1 Ф++(В) 1
(22)
;
Дадим модели (16) обозначение SSA(L)-SARIMA, модели (20) - SSA-SARIMA, а модели (22) - SSA-SARIMA + SARIMA.
Рассмотрим модели (14), (16), (20), (22) и сделаем следующие выводы. Модели Бокса-Дженкинса основываются на гипотезе, что изучаемый процесс является выходом линейного фильтра, на вход которого подан процесс белого шума, т. е. что член ряда у1 является взвешенной суммой текущего и предыдущих значений входного потока:
у, = ео + + еА-1 + е2а1-2 + . = ео + е(В) ^.
Если же мы рассмотрим модели (14), (16), (20), (22), то изучаемый процесс в таком случае будет являться уже выходом гребёнки линейных фильтров, на входы которой поданы процессы белого шума.
4. 2. 2. 2 Метод на основе совместного использования методов «Гусеница»-88А и Бокса-Дженкинса декомпозиционного подхода к прогнозированию нестационарных временных рядов
Приведём алгоритм и схему (рис. 3) метода на основе совместного использования методов «Гусени-ца»-SSA и Бокса-Дженкинса декомпозиционного подхода к прогнозированию нестационарных ВР:
1. Применяя этапы 1 и 2 метода «Гусеница»-SSA к исходному ВР определяются ВР
у 11), 1 = если используем модель (15) или (16);
у 1 , 1 = 1,г , если используем модель (19) или (20);
у 11), 1 = 1, г , е 1 = у 1 - Е у 11), если используется модель 1=1
(21) или (22)
с более простой структурой и интерпретируемые как тренды, квазипериодические и шумовые составляющие.
2. Для каждого ВР
у 11), 1 = , если используем модель (15) или (16);
у 1 , 1 = !>Т , если используем модель (19) или (20);
у!/ , 1 = 1,г , е1, если используется модель (21) или (22)
производится структурная и параметрическая идентификация моделей САРПСС (используя алгоритм раздела 4. 1. 2. 1)
У « = -
е:
у «=-
фу еу
у «=-
ф еУ
ф
(в) и -
ау , 1 = , для модели (15) или (16);
(B)
(В) в -
ау , 1 = 1,г , для модели (19) или (20);
(B)
(B) у« .г- е-(в)
(в)
ay , i = 1,r, Et =
ф ++(в)
а, для модели
(21) или (22),
используя алгоритм, приведенный в пункте 4. 1. 2. 1.
3. Вычисляются прогнозы с заданным упреждением 1 _
у 11) (1), 1 = 11Ь , для модели (15);
у 1(1), 1 = 1,г , для модели (19) или прогнозы по модели (20); _
у 1(1), 1 = 1,г , Л ( ) для модели (21) или прогнозы по модели (22).
4. Синтезируется общий прогноз как сумма частных прогнозов
L
у1 (1) = е у 1)(1), для модели (15) или вычисляется
1=1
прогноз по модели (16), предварительно приведя её к разностному выражению;
у1 (1) = е у 1)(1), для модели (19) или вычисляется
1=1
прогнозы по модели (20), предварительно приведя её к разностному выражению;
yt (1) = е у 1)(1) + ^1 (1), для модели (21) или вычисля-
1=1
ется прогнозы по модели (22), предварительно приведя её к разностному выражению.
Исходный временной ряд
Разложение в.р. методом «Гусеница*-55Д
Группировка №1 компонент разложений в р.
SARIMA-модель №1
Группировка N02 компонент
разложения в.р. SARIMA-модель
1 ill ¿1
и/ , u i .li. 1 ilU ¡At №2
Группировка Юг или SARIMA-модель
компонент разложения в.р. №1 для (15), (16); №г для (19), (20)
№(г+1) для (21),
(22)
Синтез модели (15), (16), (19), (20), (21) или (22)
Проверка адекватности синтезированной модели
Прогноз синтезированной моделью
p
p
p
Рис. 3. Схема декомпозиционного подхода к прогнозированию нестационарных временных рядов на основе совместного использования методов «Гусеница»-SSA и Бокса-Дженкинса
5. Результаты исследований
Исследование предлагаемых гибридных моделей на основе методов «Гусеница»-SSA и Бокса-Дженкинса осуществим сопоставляя их результаты прогнозов с результатами прогнозов, получаемых классическими вероятностными моделями SARIMA, обобщёнными на случай нескольких сезонных компонент. Реализация рассмотренных моделей производилась в математическом пакете MATLAB R2014a.
Тестирование будем проводить на трёх ВР потребления целевых продуктов в системах энергетики. Сначала возьмём ВР часовых данных потребления электроэнергии за трое суток (рис. 4).
Модели будем идентифицировать на выборке объёма 50 значений и производить одношаговый прогноз на 1 час вперёд, а далее повторять процедуру скольжением окна вплоть до конца ВР.
На рис. 5 и рис. 6 приведены графики автокорреляционной и частной автокорреляционной функций
1 2 3 4 5 6 7а 9 1011 12131415161718192021 222324252627232931}31 32 333435 36 37 333940414243444546474349 5051 52 53 54 55 56 57 58 596ПБ162636465666763697071727374
Рис. 4. Временной ряд часовых данных потребления электроэнергии за трое суток (в МВт)
0 2 4 6 3 1» 12 14 16 13 2» 22 24 26 23 30 32 34 36 33 40 42 44 46 43
Рис. 5. Автокорреляционная функция временного ряда потребления электроэнергии (п=50; k=48)
О 2 4 6 Б 10 12 14 16 13 20 22 24 26 28 30 32 34 26 38 40 42 44 46 43
Рис. 6. Частная автокорреляционная функция временного ряда потребления электроэнергии (п=50; k=48)
соответственно после следующего преобразования исходного ВР у, = У1у(. Анализ данных графиков совместно с отбором моделей SARIMA по критериям AIC и В1С определили следующую структуру модели:
(1 + 0.093В)(1+0.083В24) ^Ш^аОх(0,0,1)24: у( = 1 у(1 -0.476В"
Определимся с выбором параметра длины окна метода «Гусеница»-SSA. Для базового алгоритма метода «Гусе-ница»-SSA чем больше длина окна, тем более деталье ным получается разложение исходного ВР. Однако, не имеет смысла брать длину окна, большую чем половина длины ряда, т. к. сингулярные разложения траекторной матрицы, соответствующие выбору длины окна Ь и N^+1 эквивалентны. Таким образом, наиболее детальное разложение достигается при Ь=п/2 за исключением рядов конечного ранга. Маленькая длина окна может привести к смешиванию интерпретируемых компонент ряда [35]. Чтобы достичь лучшей разделимости аддитивных составляющих ВР, нужно выбирать большую длину окна.
Таким образом, выберем длину окна метода «Гусе-ница»-SSA Ь=24 и рассмотрим собственные значения траекторной матрицы (рис. 7).
Рассмотрим также графики рядов разложения (рис. 8). Анализ собственных значений траекторной матрицы, рядов разложения или факторных векторов помогает на этапе группировки так выбрать подмножества 11, 12,..., 1г, чтобы компоненты разложения были интерпретированы как трендовая, сезонная или шумовая составляющая исходного ВР.
Рис. 7. Собственные значения траекторной матрицы временного ряда потребления электроэнергии
Анализ диаграмм (рис. 9, 18, 25) говорит о снижении ошибки прогнозирования при добавлении к ЛРФ метода «Гусеница»-SSA модели SARIMA, построенной на остаточных ошибках модели, определяемых путём вычитания из значений исходного ВР соответствующих значений, получаемых при помощи ЛРФ.
На следующих двух столбчатых диаграммах (рис. 10, 11) представлены зависимости ошибок прогнозирования моделей ЛРФ-SSA + SARIMA, Тренд-88А-8АШМА и моделей ЛРФ-88А + БАШМА и ЗЗА^АШМА соответственно.
Рассмотрение приведенных выше результатов позволяет сделать вывод, что анализ собственных значений траекторной матрицы и ВР разложения исходного ВР не позволяет точно определиться с количеством отбираемых в декомпозиционную модель ВР разложения, чтобы такая модель имела наименьшую ошибку прогнозирования, а позволяет только определиться с группировкой ВР разложения для модели SSA-SARIMA. Поэтому для следующих тестовых примеров не имеет смысла приводить собственные значения траекторной матрицы и серии разложения.
Сравнительный анализ эффективности прогнозирования рассмотренными моделями будем осуществлять при помощи наиболее широко используемой статистики RMSE (Root Mean Squared Error):
RMSE = . -£(Yt - Y
где п1 - количество вычисленных прогнозов, у( - фактические значения ВР, у( - прогнозные значения ВР.
Тестирование будет проводиться следующим образом. Длина окна метода «Гусеница»-SSA будет фиксированной, но модели одного класса будут отличаться количеством первых компонент разложения, используемых в них. Поэтому на приведенных ниже столбчатых диаграммах на оси абсцисс будем откладывать количество первых компонент разложения, используемых в конкретной модели, а на оси ординат - ошибку прогнозирования, выраженную в статистике RMSE.
Анализ диаграмм (рис. 9-11, 18, 19, 25-27) позволяет сделать вывод, что увеличение количества используемых в моделях компонент разложения постепенно снижает ошибку прогнозирования, а затем, начиная с некоторого количества используемых компонент разложения, эта ошибка начинает постепенно возрастать.
Рис. 8. Графики временных рядов разложения методом «Гусеница»-5$А временного ряда потребления электроэнергии
i t=1
I
■ ЧДК1МЛ
ii
11
и
ШНШишииН
1 2 3 Л 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
Рис. 9. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA и ЛРФ-SSA + SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Рис. 10. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA + SARIMA и Тренд-SSA-SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Рис. 11. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA + SARIMA и SSA-SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Под структурной идентификацией приведенных декомпозиционных моделей будем подразумевать перебор различных вариантов разложения и выбор наилучшей по точности прогнозирования модели из множества альтернатив. Поэтому приоритетной задачей является разработка показателей, позволяющих автоматизировать и ускорить процесс определения
структуры наилучшей по точности прогнозирования декомпозиционной модели из определённого класса без перебора множества вариантов.
На следующей диаграмме (рис. 12) изображены ошибки прогнозирования RMSE для наилучших моделей из рассмотренных классов моделей. Наименьшую ошибку прогнозирования показала модель SSA-SARIMA.
Рис. 12. Диаграмма минимальных ошибок прогнозироваг ния RMSE для рассмотренных классов моделей
Как видно из диаграммы рис. 11 наилучшая по эффективности прогнозирования модель SSA-SARIMA включила в себя первые 15 компонент разложения методом «Гусеница»-SSA исходного ВР. При этом следует отметить, что группировка компонент разложения производилась как показано на следующих графиках (рис. 13).
Рис. 13. Группировка компонент разложения для наилучшей по эффективности модели SSA-SARIMA
Приведём прогноз наилучшей из класса SSA-SARIMA модели (рис. 14).
Рассмотрим более длинный ВР часовых значений потребления электроэнергии объёма 1008 значений, что соответствует 6 неделям (рис. 15). Обучение мо-
делей будет производиться на выборке данных за 5 недель (840 значений), а тестирование - на данных последней недели. Прогноз тоже будет выполняться одношаговый и производится дальнейшее скольжение окна до последнего значения ВР.
На рис. 16 и рис. 17 приведены графики автокорреляционной и частной автокорреляционной функций соответственно после следующего преобразования исходного ВР у, = У1У124у1. Анализ данных графиков совместно с отбором моделей SARIMA по критериям AIC и В1С определили следующую структуру модели: SARIMA(7,1,2) х (3,1,1)24 х (1,0,1)120 х (0,0,1)168:
Yt =
(1+0.495B + 0.947B2 )(1 - 0.861B2
v1v24 (1 + 0.396B+0.838B2 - 0.011B3 -(1 - 0.208B120 )(1 + 0.336B168) " (1 - 0.16B24 - 0.001B48 + 0.022B72 )(1 + 0.040B120)'
Выберем длину окна метода «Гусеница»-SSA Ь=396. В модели будем отбирать от 1 до 24 первых компонент разложения.
Приведенная ниже столбчатая диаграмма (рис. 18), говорит о существенном снижении ошибки прогнозирования при добавлении к ЛРФ модели SARIMA, как это было описано для модели ЛРФ-SSA+SARIMA.
Yt =
Из следующих диаграмм (рис. 19, 20) видно, что для данного тестового примера наиболее эффективный прогноз показала модель ЛРФ-SSA + SARIMA.
Приведём прогноз моделью с наименьшей ошибкой RMSE (рис. 21).
Рассмотрим теперь тестовый пример прогнозирования ВР потребления природного газа. График данного ВР среднесуточного объёма газопотребления представлен на рис. 22 (объём 758 значений). Обучение моделей будем проводить на данных за два года (объём 730 значений). Тестирование на данных за 4 недели. Прогноз будет производиться одношаговый, также как и для предыдущих двух примеров.
На рис. 23, 24 приведены графики автокорреляционной и частной автокорреляционной функций соответственно после следующего преобразования исходного ВР Y t = V1Yt. Анализ данных графиков совместно с отбором моделей SARIMA по критериям AIC и BIC определили следующую структуру модели:
SARIMA(7,1,1)x(1,0,1)21 x(1,0,0)35 x(1,0,0)56 x(1,0,0)70. (1-0.851B)
V1 (1 -0.899B+0.058B2 -0.009B3 -0.016B4+0.044B5+0.013B6 + 0.037B7) (1 - 0.609B21)
0.102B4 + 0.076B5 + 0.033B6 - 0.082B
(1-0.724B21 )(1 - 0.08B35 )(1 - 0.039B56 )(1 - 0.066B7'
Рис. 14. График временного ряда потребления электроэнергии и прогнозы, полученные моделью с наименьшей ошибкой RMSE
Рис. 15. График временного ряда часовых значений потребления электроэнергии за 6 недель (в МВт)
к— 1
иО-п^а. .л.—- Г"''1 г- " I'-- - ''-"41 "И-"' " Ч" К..,«! Г1"
г
Рис. 16. Автокорреляционная функция временного ряда потребления электроэнергии (п=840; k=200)
Рис. 17. Частная автокорреляционная функция временного ряда потребления электроэнергии (п=840; k=200)
Рис. 18. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA и ЛРФ-SSA+ SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
при добавлении к ЛРФ модели SARIMA, как это было и для двух предыдущих тестов.
Рис. 19. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA+SARIMA и SSA-SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Ниже приведена столбчатая диаграмма (рис. 25), которая говорит о снижении ошибки прогнозирования
Рис. 20. Диаграмма минимальных ошибок прогнозирования RMSE для рассмотренных классов моделей
Для случаев прогнозирования ВР со сложными стохастическими трендами модели Тренд-SSA-SARIMA дают более эффективный прогноз в сравнении с моделями ЛРФ-SSA + SARIMA (рис. 26).
Кроме того, более простая модель Тренд-SSA-SARIMA может давать более эффективный прогноз в сравнении с моделью SSA-SARIMA (рис. 27, 28).
Для ВР потребления природного газа, который имеет полиномиальный, полигармонический и стохастический тренд со сложной корреляционной структурой наилучший результат выдала модель из класса Тренд-SSA-SARIMA (рис. 28).
Ниже приведём прогнозы наилучшей из класса Тренд-SSA-SARIMA модели (рис. 29).
Таким образом, для ВР различного объёма и структуры, наиболее эффективный прогноз дают модели из различных классов.
В дальнейшем целесообразно обобщит гибридные математические модели и методы, основанные на методах «Гусеница»-SSA и Бокса-Дженкинса для прогнозирования взаимосвязанных нестационарных ВР.
Рис. 21. График временного ряда потребления электроэнергии и прогнозы, полученные наилучшей
из рассмотренных моделей
Рис. 22. Временной ряд суточных данных потребления природного газа объёмом 758 значений
(в Мм/{ут )
1 -у™ ^»"V ^Г*- 1--Г--.» 1-{,-ч Уг11 Ч^т ^ 1 Т+ гЛ**-"^-'
• А *
№ I» Ш 1« 1Н <К ¡« ¡» ¡'1 » ИО Щ ]» V) )» 1Н Щ
Рис. 23. Автокорреляционная функция временного ряда потребления природного газа (п=730; к=400)
Рис. 24. Частная автокорреляционная функция временного ряда потребления электроэнергии (п=840; к=200)
Рис. 25. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA и ЛРФ-SSA + SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Рис. 26. Диаграммы зависимости ошибки прогнозирования RMSE моделей ЛРФ-SSA + SARIMA и Тренд-SSA-SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Рис. 27. Диаграммы зависимости ошибки прогнозирования RMSE моделей Тренд-SSA-SARIMA и SSA-SARIMA от количества используемых в них первых компонент разложения метода «Гусеница»-SSA
Рис. 28. Диаграмма минимальных ошибок прогнозирования RMSE для рассмотренных классов моделей
Полученные результаты позволяют сделать заключение, что для получения эффективных прогнозов необходимо производить декомпозицию исследуемых ВР различными способами и отбирать наиболее эффективные из них для ВР определённой структуры. А также комбинировать различные модели, описывающие как статистические, так и детерминированные составляющие ВР, что позволяет достичь наилучшего качества прогнозирования. Предложенные модели учитывают следующие виды нестационарности: по тренду, по частоте и амплитуде колебаний, по дисперсии шумовой составляющей процесса.
Основным преимуществом предлагаемых методов построения моделей адекватных исследуемому процессу является их строгая формализация и, следовательно, возможность полной автоматизации всех этапов построения и использования.
Рис. 29. График потребления природного газа и прогнозы, полученные наилучшей из рассмотренных моделей
6. Выводы
П 3.
В данной статье рассмотрены трендовыи и декомпозиционный подходы к прогнозированию нестационарных ВР. В соответствии с ними предложены декомпозиционные модели различной сложности для прогнозирования нестационарных ВР, а также методы идентификации этих моделей 4.
на основе совместного использования методов «Гу-сеница^-SSA и Бокса-Дженкинса, обобщённого для случая нескольких сезонных компонент. Гибрид- 5
ные математические модели трендового подхода к прогнозированию, основанные на методах «Гу-сеница^-SSA и Бокса-Дженкинса заключаются в моделировании процесса как отклонения фактиче- 6
ских значений ВР относительно трендовой составляющей, в роли которой в предлагаемых моделях выступает линейная рекуррентная формула (ЛРФ) метода «Гусеница^-SSA или её аппроксимация мо- 7.
делью SARIMA. Основной же целью декомпозиционного подхода к прогнозированию на основе методов «Гусеница^-SSA и Бокса-Дженкинса является разложение методом «Гусеница^-SSA исходного ВР g
на множество ВР с более простой структурой, рассматриваемых независимо друг от друга; прогнозирование данных компонент разложения моделями SARIMA и вычисление общего прогноза, объединяя прогнозы построенных упрощённых моделей. 9. Предложенные модели были протестированы на ВР потребления электроэнергии и природного газа, а их результаты прогнозирования были сравнены с 10. результатами, полученными классическими вероятностными моделями SARIMA, обобщёнными на случай нескольких сезонных компонент.
Литература
1. Седов, А. В. Моделирование объектов с дискретно-распределёнными параметрами: декомпозиционный подход [Текст] / А. В. Седов. - Южный научный центр РАН. -М. : Наука, 2010. - 438 с.
2. Бэнн, Д. В. Сравнительные модели прогнозирования электрической нагрузки [Текст] / Д. В. Бэнн, Е. Д. Фармер; пер. с англ. - М. : Энергоатомиздат, 1987. - 200 с.
Qiang, Z. Singular Spectrum Analysis and ARIMA Hybrid Model for Annual Runnoff Forecasting [Text] / Z. Qiang, D. W. Ben, H. Bin, P. Yong, L. R. Ming // Water Resour Manage. - 2011. - Vol. 25, Issue 11. - P. 2683-2703. doi: 10.1007/s11269-011-9833-y
Евдокимов, А. Г. Оперативное управление потокора-спределением в инженерных сетях [Текст] / А. Г. Евдокимов, А. Д. Тевяшев. - Х.: Вища школа, 1980. - 144 с. Lawrance, A. J. Stochasting modelling of overflow time series [Text] / A. J. Lawrance, N. T. Kottegoda // J. R. Stat. Soc. A. - 1977. - Vol. 140, issue 1. - P. 1-47. doi: 10.2307/2344516 Fernando, D. A. K. Generation and forecasting of monsoon rainfall data [Text] / D. A. K. Fernando, W. A. Jayawardena // In Proc. of the 20th WEDC conference. Colombo, Sri Lanka, 1994. - P. 310-313.
Yurekli, K. Application of linear stochastic models to monthly flow data of Kelkit Stream [Text] / K. Yurekli, A. Kurunca, F. Ozturkb // Ecol Model. - 2005. - Vol. 183, Issue 1. - P. 67-75. doi: 10.1016/j.ecolmodel.2004.08.001 Broomhead, D. S. Extracting qualitative dynamics from experimental data [Text] / D. S. Broomhead, G. P. King // Physica D. - 1986. - Vol. 20, Issue 2-3. - P. 217-236. doi: 10.1016/0167-2789(86)90031-x Fraedrich, K. Estimating the dimension of weather and climate attractor [Text] / K. Fraedrich // J. Atmos Sci. -1986. - Vol. 43. - P. 419-432.
Vautard, R. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series [Text] / R. Vautard, M. Ghil // Physica D. - 1989. - Vol. 35, Issue 3. - P. 395-424. doi: 10.1016/0167-2789(89)90077-8
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
i
Ghil, M. Interdecadal oscillations and the warming trend in global temperature time series [Text] / M. Ghil, R. Vautard // Nature. -1991. - Vol. 350, Issue 6316. - P. 324-327. doi: 10.1038/350324a0
Yiou, P. Spectral analysis of climate data [Text] / P. Yiou, E. Baert, M.F. Loutre // Surv Geophys. - 1996. - Vol. 17, Issue 6. -P. 619-663. doi: 10.1007/bf01931784
Lisi, F. Combination of singular spectrum analysis and auto regressive model for short term load forecasting [Text] / F. Lisi, O. Nicolis, M. Sandri // Neural Process Lett. - 1995. - Vol. 2, Issue 4. - P. 6-10.
Sivapragasam, C. Rainfall and discharge forecasting with SSA-SVM approach [Text] / C. Sivapragasam, S.Y. Liong, M.F.K. Pasha // J. Hydroinform. - 2001. - Vol. 3, Issue 7. - P. 141-152.
Golyandina, N. Analysis of time series structure: SSA and related techniques [Text] / N. Golyandina, V. Nekrutkin, A. Zhigljavsky // Chapman and Hall/CRC, New York, 2001. doi: 10.1201/9781420035841
Marques, C. A. F. Singular spectrum analysis and forecasting of hydrological time series [Text] / C. A. F. Marques, J. A. Ferreira, A. Rocha, J. M. Castanheira // In Meeting of the European-Union-of-Geosciences. Vienna, Austria, 2005.
Hassani, H. Forecasting European industrial production with singular spectrum analysis [Text] / H. Hassani, S. Heravi,
A. Zhigljavscky // Int. J. Forecast. - 2009. - Vol. 25, Issue 1. - P. 103-118. doi: 10.1016/j.ijforecast.2008.09.007
Wensheng, D. Financial Time Series Forecasting Using A Compound Model Based on Wavelet Frame and Support Vector Regression [Text] / D. Wensheng, Lu Chi-Jie // In the 4th International Conference on Neural Computation, 2008. - P. 328-332. doi: 10.1109/icnc.2008.455
Kurbatskii, V. G. On the Neural Network Approach for Forecasting of Nonstationary Time Series on The Basis of the Hilbert-Huang Transform [Text] / V. G. Kurbatskii, D. N. Sidorov, V. A. Spiryaev, N. V. Tomin // Automation and Remote Control. - 2011. -Vol. 72, Issue 7. - P. 1405-1414. doi: 10.1134/s0005117911070083
Zhang, W. Q. Time series forecasting method based on Huang transform and BP neural network [Text] / W. Q. Zhang, C. Xu // In Proc. Of the 7th International Conference on Computational Intelligence and Security, 2011. - P. 497-502. doi: 10.1109/cis.2011.116 Lu, C.-J. ICA-Based Signal Reconstruction Scheme with Neural Network in Time Series Forecasting [Text] / C.-J. Lu, J.-Yu Wu, T.-S. Lee // In First Conference on Intelligent Information and Database Systems, 2009. - P. 318-323. doi: 10.1109/aciids.2009.28 Xiang, L. A hybrid support vector regression for time series forecasting [Text] / L. Xiang, Y. Zhu, G.-J. Tang // In World Congress on Software Engineering, 2009. - P. 161-165. doi: 10.1109/wcse.2009.130
Sallehuddin, R. Hybridization Model of Linear and Nonlinear Time Series Data for Forecasting [Text] / R. Sallehuddin, S. M. Shamsuddin, S. Z. M. Hashim // In Second Asia International Conference on Modelling & Simulation, 2008. - P. 597-602. doi: 10.1109/ams.2008.142 Henrique, S. H. Combining Neural Networks and ARIMA Models for Hourly Temperature Forecast [Text] / H. S. Hippert, C. E. Pedreira, R. C. Souza // In Proceedings of the International Joint Conference on Neural Networks, 2000. - P. 1-6. doi: 10.1109/ijcnn.2000.860807 Xuemei, L. A Novel Air-conditioning Load Prediction Based on ARIMA and BPNN Model [Text] / L. Xuemei, D. Lixing, S. Ming, X. Gang, L. Jibin // In Asia-Pacific Conference on Information Processing, 2009. - P. 51-54. doi: 10.1109/apcip.2009.21 Tian, F. P. Forecast of Cerebral Infraction Incidence Rate Based on BP Neural Network and ARIMA Combined Model [Text] / F. P. Tian, L. L. Ma //In International Symposium on Intelligence Information Processing and Trusted Computing, 2010. -P. 216-219. doi: 10.1109/iptc.2010.7
Feng, K. Time Series Forecasting Model with Error Correction by Structure Adaptive Support Vector Machine [Text] / K. Feng, W. Xiaojuan // In International Conference on Computer Science and Software Engineering, 2008. - P. 1067-1070. doi: 10.1109/csse.2008.88 Lo, J.-H. A Data-Driven Model for Software Reliability Prediction [Text] / J.-H. Lo // In International Conference on Granular Computing, 2012. - P. 1-6. doi: 10.1109/grc.2012.6468581
He, Y. Research on Hybrid ARIMA and Support Vector Machine Model in Short Term Load Forecasting [Text] / Y. He, Y. Zhu, D. Duan // In Proceedings of the Sixth International Conference on Intelligent Systems Design and Applications, 2006. - P. 1-5. doi: 10.1109/isda.2006.229
Linh, B. N. An Empirical Study on Forecasting using Decomposed Arrival Data of an Enterprise Computing System [Text] /
B. N. Linh, A. Amy, H. Doug // In 9th International Conference on Information Technology- New Generations, 2012. - P. 756-763. doi: 10.1109/itng.2012.36
Hou, Z. Standardized Software for Wind Load Forecast Error Analyses and Predictions Based on Wavelet-ARIMA Models -Applications at Multiple Geographically Distributed Wind Farms [Text] / Z. Hou, Y. V. Makarov, N. A. Samaan, P. V. Etingov // In Hawaii International Conference on System Sciences, 2013. - P. 5005-5011. doi: 10.1109/hicss.2013.495
Щелкалин, В. Н. Трендовый подход прогнозирования временных рядов на основе метода «Гусеница^-SSA [Текст] / Материалы 14-й Международной научно-технической конференции SAIT 2012, Киев, 24 апреля 2012 г. / В .Н. Щелкалин // УНК "ИПСА" НТУУ "КПИ". - К.: УНК "ИПСА" НТУУ "КПИ", 2012. - С. 258 - 259.
Vahabie, A. H. Combination of Singular Spectrum Analysis and Autoregressive Model for short term load forecasting [Text] / A. H. Vahabie, M. M. R. Yousefi, B. N. Araabi, C. Lucas, S. Barghinia // IEEE LAUSANNE POWERTECH, 2007. - P. 1090-1093. doi: 10.1109/pct.2007.4538467
Щелкалин, В. Н. Декомпозиционный подход прогнозирования временных рядов на основе метода «Гусеница^-SSA [Текст] : матер. 14-й Междун. науч.-тех. конф. SAIT / В. Н. Щелкалин // УНК "ИПСА" НТУУ "КПИ". - К.: УНК "ИПСА" НТУУ "КПИ", 2012. - С. 260-261.
Голяндина, Н. Э. Метод «Гусеница^-SSA: прогноз временных рядов [Текст]: уч. пос. / Н. Э. Голяндина. - СПб. : С.-Петербургский государственный университет, 2004. - 52 с.
