Гибридные математические модели и методы прогнозирования временных рядов с учётом внешних факторов Текст научной статьи по специальности «Математика»

-□ □-

В роботi запропоновано математичн моделi прогно-зування нестащонарних часовихряЫв зурахуванням зов-тштх чиннитв i методи гх структурног идентифжацп, засноват на стльному використанн багатовимiрного варiанта методу «Гусениця»-ББЛ та моделей БЛЯ1МЛХ, розширених для прогнозування часових рядiв з кшько-ма перюдичними компонентами i з урахуванням декшь-кох екзогенних змтних. Експериментальн результати показують високу ефективтсть запропонованих моделей прогнозування при виборi видповидних структурних параметрiв у порiвнянш з моделями БЛЯМЛХ

Ключовi слова: прогнозування, структурна и)ен-тифЫащя, декомпозицшна модель, метод Бокса-

Дженктса, метод «Гусениця»-ББА

□-□

В работе предложены математические модели прогнозирования нестационарных временных рядов с учётом внешних факторов и методы их структурной идентификации, основанные на совместном использовании многомерного варианта метода «Гусеница»-ББА и моделей БАЯМАХ, расширенных для прогнозирования временных рядов с несколькими сезонными составляющими и с учётом нескольких экзогенных переменных. Экспериментальные результаты показывают высокую эффективность предложенных моделей прогнозирования при выборе подходящих структурных параметров в сравнении с моделями БАШМАХ

Ключевые слова: прогнозирование, структурная идентификация, декомпозиционная модель, метод Бокса-

Дженкинса, метод «Гусеница»-ББА -□ □-

УДК 519.254

|DOI: 10.15587/1729-4061.2014.31729|

ГИБРИДНЫЕ МАТЕМАТИЧЕСКИЕ МОДЕЛИ И МЕТОДЫ ПРОГНОЗИРОВАНИЯ ВРЕМЕННЫХ РЯДОВ С УЧЁТОМ ВНЕШНИХ ФАКТОРОВ

В. Н. Щелкалин

Инженер

Кафедра прикладной математики Харьковский национальный университет радиоэлектроники пр. Ленина, 14, г. Харьков, Украина, 61166 E-mail: vitalii.shchelkalin@gmail.com

1. Введение

Точность решения задач оперативного планирования режимов распределения целевых продуктов (воды, газа, электроэнергии) в инженерных сетях, а также возможность использования полученных результатов в практике диспетчерских служб управления инженерными сетями, определяются двумя факторами: степенью адекватности используемых математических моделей и точностью исходных данных. Важнейшие исходные данные для задач оперативного планирования - это прогнозы потребления целевого продукта в инженерных сетях. Недостаточное знание механизма процессов водо-, газо- и электропотребления, влияние на них огромного числа неконтролируемых и неуправляемых факторов приводят к тому, что эти процессы являются случайными процессами, которые в зависимости от типа потребителя могут носить как стационарный, так и нестационарный характер. Кроме того, они подвержены влиянию метеорологических, организационных и хронологических факторов, изменяющих их стохастическую структуру [1]. Это привело к появлению большого числа методов прогнозирования этих процессов, большинство из которых носит частный характер и позволяет получить удовлетворительные результаты только для ограниченного класса потребителей и типов прогнозов.

В настоящей работе в качестве моделей процессов потребления целевых продуктов (ППЦП) в инженерных сетях рассмотрим гибридные модели на основе многомерного варианта метода «Гусеница»-$$А (MSSA) и сезонных моделей авторегрессии - проинтегрированного скользящего среднего с экзогенными переменными, расширенных на случай прогнозирования ВР с несколькими сезонными компонентами и с учётом нескольких экзогенных переменных (далее такие расширенные модели будем обозначать САРПС-СЭ, на англ. - SARIMAX), позволяющих получить адекватные модели ППЦП, реальных потребителей с различной стохастической, детерминированной или вероятностно-детерминированной структурой.

Как уже отмечалось, ППЦП в инженерных сетях зависят от большого числа факторов трех основных типов: хронологических (время года, дни недели, сменность работы, режимы работы промышленных предприятий и т. п.), метеорологических (температура окружающей среды, скорость ветра, влажность наружного воздуха и т. п.), организационных (мероприятия, проводимые в сетевых системах, массовые мероприятия, проводимые с населением).

Соотношение между компонентами ППЦП определяется типом потребителя. Влияние хронологических факторов на ППЦП приводит, как правило, к появлению в последних периодических компонент. Их природа связана с характером технологических

©

процессов промышленных предприятии, временами года и т. п. Таким образом, появляется возможность для построения широкого класса моделеИ ППЦП, содержащих периодические компоненты либо в виде одноИ синусоиды, либо в виде отрезка ряда Фурье, либо в виде амплитудно-модулированного процесса. Влияние хронологических факторов для некоторых типов потребителей, а точнее, совместное влияние всех факторов может уменьшить или даже полностью исключить периодические компоненты из исследуемых процессов. Непрерывное увеличение потребления воды, газа или электроэнергии, а также влияние сезонных компонент с большими периодами приводит к появлению в этих процессах детерминированных и стохастических трендов и, следовательно, к нестационарности ППЦП. Установление этого факта позволило использовать модели и методы, позволяющие описывать ППЦП как нестационарные процессы.

Влияние на ППЦП температуры окружающей среды послужило основой для построения моделеИ, описывающих процессы газо- и электропотребления как некоторую детерминированную функцию от температуры окружающей среды.

В методе Бокса-Дженкинса для связи параметров ВР с экзогенными факторами используется мультипликативная модель SARIMAX, обладающая рациональной структурой и позволяющая адекватно описывать как стационарные, так и однородные нестационарные случайные процессы, содержащие nS периодических компонент. Для связи параметров ВР используются модели линейных дискретных передаточных функций с запаздывающим аргументом. Рассмотрены методы проверки адекватности получаемых моделей, основанные на анализе автокорреляционных и взаимных корреляционных функций остаточных ошибок модели. Приведены алгоритмы вычисления и коррекции прогнозов и их дисперсий.

Отличительной чертой метода MSSA является то, что он не требует предварительного задания модели ряда. Тем не менее, он позволяет раскладывать временной ряд (ВР) на интерпретируемые составляющие, такие как тренд (в частности, линейный или экспоненциальный), периодические компоненты и шум. При этом не надо заранее знать параметрический вид тренда, а также о наличие колебательных компонент и их периодах [2]. Техническую основу метода составляет сингулярное разложение траекторной матрицы, столбцами которой являются вектора вложения - отрезки исходных ВР длины L, основного параметра метода, называемого длиной окна [3]. Анализ членов сингулярного разложения позволяет сначала классифицировать их как относящиеся к одной из компонент ряда, а затем выделить эту компоненту.

ВР потребления целевых продуктов в инженерных сетях являются существенно нелинейными со средне- и долгосрочными зависимостями, поэтому аппроксимация таких сложных ВР только моделью АШМА не всегда является удовлетворительной. Сингулярный спектральный анализ (SSA) является одним из наиболее эффективных методов для анализа ВР со сложными периодическими компонентами. SSA стал широко используемым методом анализа периодических и трендовых составляющих ВР [8-11]. В последнее время были разработаны новые модели для прогнозирования, основанные на SSA. В [12] предложили линейную рекуррентную формулу сингулярного спектрального анализа (SSA-LRF), и эта модель была использована для решения некоторых практических задач [13, 14].

В [15] предложена комбинированная аддитивная модель прогнозирования, в которой к ЛРФ добавляется модель авторегрессии - проинтегрированного скользящего среднего. Исследования, комбинирующие метод «Гусеница»-SSA с другими моделями можно встретить и в зарубежной литературе. Так, например, в [16] разработана подобная модель, в которой ЛРФ комбинируется с моделью авторегрессии и применяется для краткосрочного прогнозирования нагрузки иранской национальной энергосистемы. Другой пример совместного использования метода «Гусеница»-SSA с моделями Бокса-Дженкинса можно найти в [17]. Там методом «Гусеница-SSA производится разложение исходного ВР, затем производится группировка ВР разложения на интерпретируемые составляющие, такие как трендовая, сезонные и шумовые. Далее для каждого сгруппированного ВР идентифицируется и оценивается модель ARIMA. Остаточная ошибка моделирования корректируется авторегрессионной моделью и определяются прогнозы синтезированной моделью. Примерно в то же время в [18] была предложена подобная модель, однако, в отличие от [17], для идентификации и оценивания сгруппированных ВР использовались сезонные модели авторегрессии - проинтегрированного скользящего среднего (САРПСС, на англ. - SARIMA), а также такая гибридная модель была обобщена уже для случая прогнозирования ВР, зависящего от нескольких экзогенных ВР. Гибридные модели на основе метода SSA, моделей ARIMA и искусственных нейронных сетей активно обсуждаются в [42-45].

Таким образом, анализ рассмотренной литературы позволяет сделать вывод, что в настоящее время происходит отход от статистической постановки задачи прогнозирования ВР и эффективное прогнозирование нестационарных ВР предполагает использование различных приёмов декомпозиции, синтеза прогнозных моделей и отбора их из множества альтернатив.

2. Анализ литературных данных и постановка проблемы

3. Цели и задачи исследования

В последнее время появилось много работ, посвященных гибридным декомпозиционным моделям и методам прогнозирования, краткий обзор которых представлен в [4]. Для оценивания эффективности прогнозы обычно сравниваются с прогнозами, получаемыми моделями ARIMA [5-7].

Целью проведенных исследований является разработка вероятностно-детерминированных и декомпозиционных математических моделей и методов их идентификации для прогнозирования нестационарных ВР со сложной структурой с учётом экзогенных факторов и проверка эффективности прогнозирования

разработанными моделями в сравнении с прогнозами, получаемыми вероятностными моделями SARIMAX, обобщёнными для прогнозирования ВР с несколькими сезонными составляющими и несколькими экзогенными переменными [1].

Для этих целей в работе использованы метод MSSA и модели SARIMAX. На их основе разработаны гибридные математические модели и методы прогнозирования нестационарных ВР с несколькими сезонными составляющими с учётом нескольких экзогенных факторов. Метод MSSA в работе применяется для разложения исходных прогнозируемого и экзогенных ВР на ВР с более простой структурой, а модели SARIMAX - для построения математических моделей данных компонент разложения и вычисления их прогнозов. Выбор различных параметров метода MSSA приводит к созданию, в сочетании с моделями SARIMAX, различных вариантов методов прогнозирования и к синтезу на их основе гибридных математических моделей с различными структурами.

4. Математические модели и методы прогнозирования нестационарных временных рядов с экзогенными переменными

В данном разделе сначала приведём краткое изложение методов и моделей прогнозирования нестационарных ВР, учитывающих внешние факторы: многомерного варианта метода «Гусеница»-SSA и модели SARIMAX, расширенной для прогнозирования ВР с учётом нескольких экзогенных переменных. Приведенные методы и модели необходимы далее для описания гибридных моделей на их основе.

4. 1. Модель авторегрессии - проинтегрированного скользящего среднего с экзогенными переменными

В [1] была предложена математическая модель многосвязных нестационарных случайных процессов, позволяющая описывать широкий класс регулярных случайных процессов, содержащих полиномиальные, полигармонические и стохастические тренды и аддитивный цветной шум:

у =ь<(£)хо +^Б)а у ь 8(;)(б) х'-ь'+ф++(в)

(1)

где , 1 = 1,п - исходный или преобразованный (нормированный или прологарифмированный) ВР; п - объём выборки; х[;), 1 = 1,п, ; = 1,N - исходные или преобразованные (нормированные или прологарифмированные) ВР изменения значений внешних (экзогенных) факторов; N - количество экзогенных факторов; Б - оператор сдвига по времени на одну единицу назад, такой что = хы ; Ь; - величина задержки ;-го экзогенного ВР по времени относительно прогнозируемого ВР ; §[;) (Б) и юС'(Б) - полиномы передаточных функций от Б степеней г; и с; соответственно; 8<;)(Б) = 1 -§1;)Б-8<;)Б2-...-8(;)Бг; ; ю<;)(Б) = ю0;)-ю1;)Б--ю^-Б"2 -...-ю[;)Бс;; Ф++ (Б) - обобщенныйоператор авторегрессии (АР) порядка р+= р* + Ь D;S; , р* = Ь р;Д; ;

;=1 ;=1

Ф+++(Б) = Фр*(Б)У!?1 V!?2 ; D;, ; = - порядок взя-

тия разности ; = 1,п5 - период ;-й периодической

компоненты, причем Д1=1; п5 - количество периодических компонент; VS и БД; - упрощающие операторы такие, что VSyt = (1 -БД;)у1 = у1 -у^ ; Ф** (Б) - обобщенный оператор АР порядка р* вида Фр*(Б) = ПфР;'(БД;);

ф'р'(бД; ), ; = 1,п5 - полиномы от БД; степеней р; соответственно, определяющие составляющие АР периодических компонент с периодами соответственно; 8р, (Б) - обобщенный оператор скользящего среднего

П п,

(СС) порядка я* = Ьвида 6*. (В) = Пб^В8,)

;=1 ( 1=0

б^Б^), ; = 1,п5 - полиномы от БД; степеней соответственно, определяющие составляющие СС периодических компонент с периодами соответственно; а1 - остаточные ошибки модели (1).

4. 1. 1. Проверка адекватности модели авторегрессии - проинтегрированного скользящего среднего с экзогенными переменными

Прежде чем воспользоваться моделью (1), необходимо проверить, является ли эта модель адекватной исследуемому процессу. Для этого необходимо установить, значимо ли отличаются остаточные ошибки модели от процесса "белого шума", значима ли величина взаимной корреляции между, приведенными к стационарному виду при помощи (26), процессами х^' и а1. Для принятия решения об адекватности модели используем метод, предложенный в работе [23].

Если структура модели (1) принята правильно и в соответствующие уравнения подставлены истинные значения параметров, остаточные ошибки будут "белым шумом", а выборочные корреляции будут распределены взаимонезависимо с нулевым средним зна-

- 1 , *

чением и дисперсией, равной — , где М = п - и - р ,

N N. ,'.

Если вме-

где п' = п-Ь, и = 1

Ь Г, Ь(С; + Ь; )

V ;=1 ;=1

сто истинных значений параметров подставлены их выборочные оценки, распределение корреляций с малыми задержками искажается. В частности, дисперсия выборочных корреляций с малыми задержками

может быть значительно меньше —, а их значения

М'

сильно коррелированы. Поэтому простое сравнение выборочных автокорреляций гаа(к) при малых к со

"стандартной ошибкой"

1

М'

может создать ложное

впечатление о малости первых. Кроме того, всплески выборочной автокорреляционной функции при малых задержках могут происходить просто из-за высокой корреляции между выборочными оценками. Если ам-

плитуда таких всплесков мала по сравнению

1

, они

могут быть чисто случайными и не свидетельствовать о каком-либо реальном поведении теоретических автокорреляций. Все это говорит о том, что предпочтительнее не рассматривать отдельные гаа(к), а оценить

неадекватность модели по совокупности первых автокорреляций. Можно показать [23], что если К достаточно велико, то для адекватной модели величина

е = * I ^(к)

ной функции

е;,(в)

с (в) С

и не зависит от —т^-—-. Сравнив

S = М'1 Г,2й(к)

ю'

'/(В)

8'(В)у1 = Хю(1)* (В) х^, + ф(В)а<,

где

8'(В) = И8Г!) (В) Ф++ (Вр+),

ю(')*(В) = юС')(В) Ф++ (Вр+), ; = 17ы,

Ф(В)=П ««(В) ед,(вд").

(5)

(2)

распределена примерно как %2 с К-р* -д* степенями свободы, т. е. зависит только от параметров передаточ-

Предположим, что адекватная статистическая модель процесса х['), ; = 1,К может быть экономично представлена в классе общей мультипликативной сезонной модели авторегрессии - проинтегрированного скользящего среднего, т. е.

х(')=_

/т-Л 11 OrtI5.HU.Hl 1Л /Л . .

Ф++(В) 8г;)(в)

значение Ц с таблицей значений %2 при заданном уровне значимости, можно приближенно проверить гипотезу об адекватности полученной модели.

Для проверки адекватности модели передаточной

Ю?(В)

функции д(;)(в) используем статистику, основанную

на взаимной корреляционной функции г й (к) преобразованного процесса а[;) = Ф(В)хи остаточных ошибок модели, где Ф(В) - оператор метода предварительного выравнивания спектра, который позволяет преобразовать коррелированный процесс х[;) в некоррелированный - белый шум а[;). Величина

е,(;) (В4 )

Ф+(;) (Вр+)

а[;), ; = 1,К,

(6)

где а[;), ; = 1,К - процесс типа белого шума. Подставляя (6) в (1) и обозначая

Ь(0(В) = -

;(В)еХй(ВЧ*]В

¥(В) =

8Г°(В)Ф+Й (Вр+

р+

е',(Вч")

д

Ф++ (Вр+)'

(7)

(8)

представим модель (1) в виде

(3)

у< = 1 Х(;)(В)а« + ¥(В)а< = ЦX« а« + 1^ аи ,

приближенно распределена так же, как %2 с

N N

К +1 -1(г; + с;) +1 степенями свободы, где 1(г + с; ) + 1 -

число параметров в подгоняемой модели передаточ-

ной функции (А , . . Принятие решения об адекватно-

ЬГ/ (В)

сти модели осуществляется путем сравнения значения S со значением %2 при заданном уровне значимости.

В автоматизированных системах вид выборочных функций г и (к) и гаа(к) в ряде случаев позволяет ЛПР определить причину неадекватности моделей

юС0(В) ед.(В) „

, ф+ (в) и обеспечить целенаправленный поиск

адекватных моделей. В автоматических системах поиск адекватной модели может быть выполнен в ходе решения оптимизационной задачи поиска тт^) и тт(О) в ограниченном классе рациональных структур модели.

4. 1. 2. Вычисление прогнозов временных рядов моделями сезонной авторегрессии - проинтегрированного скользящего среднего с экзогенными переменными

Запишем математическую модель (1) в более удобной форме:

=1 ]=<

где а[;) и а1 статистически независимы V ; = 1^.

Прогноз у,.(1) величины у^ будем искать в виде суммы линейных функций вида

у1(1) =1 Iа«+1^аи.

¡=1 ^о j=o

Дисперсия ошибки прогноза с упреждением I равна

Е[уь+1 - у^!)]2 =

N (г /л"|2

1([^]2+[*«]2+...+[*« ]2

+(1+¥?+...+^2-1)о2+ц оав+ ¡=1 j=0

+1 (¥н-¥н)2 ®2.

(9)

(4)

Минимум выражения (9) достигается, если Хщ = = Таким образом, прогноз у,.(1)

величины уы в момент 1 с минимальной среднеква-дратической ошибкой определяется как условное математическое ожидание уы в момент 1.

Обозначая, как и ранее, условные математические ожидания в момент времени 1 квадратными скобками, непосредственно из (4) получаем выражение для прогноза с упреждением 1:

N

р++Ь

у^1) = I ]-)=1

11' ю«"[х11ь;-,]- I Ф1 ] ,

;=1 д=0, ^ "" ^=1

где 8*, ю);» , ф) определяются из уравнения (5), а

[ )J кО), ) > °;

|х(1+),) < °; . —

-4+)> °, )> °.

Прогнозы х('»() получают в соответствии с уравнением (6). Значения а(+), ) < ° могут быть взяты из (1) в три этапа. к

1. Определяем значения у( для (>Ь(Ь; + С) +1

;=1

у (= Ь Й^Ь (8^ + 8^, + ...+ 8^,)+

;=1 8; (Б) ;=1

N ,

+Ь (ш°М-Ь;-Ю(;)х(-Ь; -1 -...-ш«х(-Ь;-с;).

Предшествующие значения у( приравниваем к нулю.

2. Находим значения

е (= у(- у (.

3. Непосредственно вычисление а( = е(п) осуществляется итерационно по формуле

8р» (Б) =8Г;'(Б)Ф+и (Бх1;)) =

Чй

= 1 -81; '**Б -82;)**Б2 -... -8р'**БР',

где Pi= Рх(;)- г; .

Для этого используем формулу

8**=-Ь 8(Ч> ,

(14)

где 8°) = -1 , Ф+(;, =-1 , 8(;) = °|) > г; , Ф+,;, = > р^о .

Вывод формул (12) и (14) достаточно прост и поэтому не приводится. Подставляя (11) и (13) в (7), имеем

(1 -81' '**Б -82; '**Б2 - . -8р'**БР' )(1 + А)' 'Б + Х 'Б2 + .) =

= (ю°'**БЬ' -ю<;'**БЬ'+1 -ю2;'**БЬ'+2 - . -ю^'). (15)

Если Ь; =°, то из (15) получаем т(; '**, ) = °; Ь8(; '**)-ю); '**,) >1,

где 8)''** = > р', ю);'** = > я' V ; = 1,К. Если Ь; >°, то

А<; Ч

Л )**

) < Ьг ) = Ь;

^ 8 '** ХМ -Ю(-Ь*,) >Ь;,

(16)

(17)

к=1 )=1

где е(°' = е( -ц , ( >Ь (Р; + ?;)Д;, ц - среднее значение ряда е(.

При вычислении дисперсии прогнозов с упреждением I необходимо определить веса А,)'» и у ■. Сначала

); )** );)** ); )** 3

находим параметры ю° , ю1' , ..., юя; оператора ю«**(Б) = ю£>(В)6*й (Бч*х);))БЬ; =

= ю°'**БЬ; -ю(;'**БЬ;+1 -ю2'гБЬ'+2-...-ю^ БЬ'+4-', (11) где я' = я*н + с; . Для этого используем формулу

ю);)** =

ю°;), 1 = °; у ю);»6*);), 1 < 1 < с;;

С;

Ью)Чв , с; <1 < Ч'.

(12)

Далее определяем параметры 8);» , 82'' , ..., 8'Р,' опе-

ратора

где 8);'** = > р', ю);'** = °|) > я'.

Вес оператора (8) находим по формуле

У) = 1

1,

) = 1;

ЬФ+У)-; ^ >1,

(18)

где Ф + = ° , ; > р*; 6* = ° , ) > я*.

4. 1. 3. Вычисление дисперсии прогнозов

Значение весов А);» и у) позволяет найти дисперсию прогноза с упреждением 1. При вычислении прогнозов по уравнению (1°) дисперсия прогнозов (9) с упреждением 1 равна

V«=Ь Ь[А);)]<й+С2Ьу2,

(19)

=1 )=1

где веса А);» и у) могут быть получены в соответствии

с уравнениями (16)-(18). Значения оценок дисперсии

22

о );) и оа получают в результате решения е метрической оптимизации моделей (6), (1).

4. 1. 4. Коррекция прогнозов

Знание весов А);» и У) позволяет не только вычислять дисперсию прогнозов, но и произвести их коррек-

цию внутри интервала упреждения Цу. При вводе в систему новых значений реализаций процессов у1+1 и х[+4 можно вычислить фактические ошибки прогнозов и осуществить их коррекцию в момент времени 1+1 для упреждений 1=1, 2, ..., Цу -1.

Коррекция прогнозов, вычисляемых в соответствии с уравнением (10), несколько сложнее. Для получения формулы коррекции прогнозов представим уравнение (10) в виде

У

(1)=х х«+х

=1 ]=|

В этом случае

Уь+1(1)=X (^1,)«[+1++^«Е+• ■ ■) ■+

N

(1 + 1) =Х (^(+)1«1') + ^(+)2«1°1 + ^2 +■■)-

Вычитая (21) из (20), имеем

(1) = У4 (1 + 1) + Х^|,)«1,) + ¥а1+1.

(20)

(21)

(22)

Это уравнение означает, что прогноз величины у1+1+1, вычисляемый в момент времени 1+1, равен прогнозу этой величины в момент времени 1 плюс ошибки прогноза на шаг вперед и а1+1 с коэффициентами

Уравнение (22) в дальнейшем используем для коррекции прогнозов. Величина а[+1 по определению равна

а 1+1 = х 1+1 - У 1,)(1).

а1+1 находим с помощью уравнения а1+1 = у 1+1- У Л1^

где

(23)

(24)

У1+1 = У1-

(,)

1=1 8Г') (В)

Из предположения, что ошибки а(,) и а подчиняются нормальному закону распределения, следует, что при известных значениях процессов х[') и уь до момента времени t условное распределение вероятности Р(У1+1| У1,У1-1,—,х1'),х['-1) будущего значения процесса У1+1 будет также нормальным с математическим ожиданием Ук(1) и дисперсией V(1) . Вероятностные пределы уы(1) для величины у^ имеют вид

уы(1) = У1(1) ± ие/ ^Ту^.

(25)

где ие/2 квантиль уровня 1 стандартного нормального распределения. При вычислении прогнозов по (10) значения У(1) в (25) вычисляются соответственно по (19). Численные значения ие/2 выбираются

равными 0.68, 1.65, 1.96 или 2.58 в зависимости от того, находятся ли будущие значения у^ в интервале У,(1)±ие/2>/У(1) с вероятностью 0.5, 0.9, 0.95 или 0.99 соответственно.

4. 1. 5. Алгоритм построения модели Бокса-Джен-кинса, расширенной для прогнозирования временных рядов с несколькими сезонными составляющими и с учётом нескольких экзогенных переменных

Данный алгоритм основан на использовании в качестве математической модели ВР выражения (1). Блок-схема алгоритма представлена на рис. 1, 2.

Алгоритм состоит из 29 этапов:

1) начало;

2) ввод обучающих выборок у^ х['), 1 = 1,^ 1 = 1,п;

3) идентификация периодов S1, S2, ..., S и порядка взятий разностей D1, Э2, ..., Dn . На основании анализа вида автокорреляционной функции процесса определяем порядок взятия разности с периодом Sn =1. Используя алгоритм выявления скрытой периодичности [24], определяем периоды периодических компонент S1, S2, ..., Sn -1. На основании анализа автокорреляционных функций процессов

У'1 = у'1 = ^^ , ..., у'1 = у5с; х„. хУ^

определяем порядки взятия разностей Di , 1,п5 -1;

4) преобразование рядов у, и х['), , = 1,N. По известным значениям Di, Si , , = 1,п5 осуществляем преобразование исходных выборок у1, х['), , = 1,N в

У'1 = VSD22 х„.хУ^,

х1')'^х„.хVD:sx(');

(26)

5) построена ли модель ряда х'1? Если «да», то переходим к 10-му этапу, если «нет» - к 6-му этапу.

Этапы 6-9-й соответствуют этапам 4-7-му алгоритма, приведеного в [4] и обеспечивают построение адекватной модели процесса х[')', , = 1^. Идентификацию модели процесса х[')', , = осуществляем в классе общей мультипликативной модели авторегрессии - проинтегрированного скользящего среднего, однако, можно использовать и другие модели, представленные в [4];

10) надо ли вычислять прогнозы х[')(1) , , = 1^? Если «да», то переходим к 11-му этапу, если «нет» - к 12-му этапу;

11) модель х['), , = 1,N адекватна? Если «да», то переходим к 13-му этапу, если «нет» - к 8-му этапу;

12) вычисление прогнозов х[')(1), 1 = 1,ЦхП , , = 1^. Осуществляем в соответствии с разностным уравнением, приведенным в [1]; --

13) ввод прогнозов х 1')(1), 1 = 1,Цх' , , =

14) модель У1 адекватна? Если «да», то переходим к 19-му этапу, если «нет» - к 15-му этапу; ю(,) (В)

15) идентификация параметров моделей (А , . ,

- , ^ «>№ МВ , = 1,N ( У (В) = Х ) ). Производится методом выравнивания спектра компоненты процесса у[,), , = 1,N в соответствии с видом и значениями параметров

модели процесса x^, i = 1,N . Идентификация параме-• ф V (в)

тров передаточной функции V3 (В) =

ф ++(В)

выполня-

ется в соответствии с алгоритмом, приведенным в [4];

16) грубая оценка параметров передаточной функции. Осуществляется путём решения системы линейных уравнений, представленных в [1]. В качестве грубых оценок параметров модели передаточной функции

в;.(в)

ф+ (в) используем константы типа ±0,1;

ошибок моделей а1 и а[1), 1 = . Проводится путём отбрасывания их Ly и Lx , 1 = соответственно старых значений и прибавлением Ly и Lx , 1 = новых значений этих рядов;

28) модель у1 адекватна? Если «да», то переходим к 10-му этапу, если «нет», то, предполагая, что новые данные привели к изменению только параметров модели, переходим к 17-му этапу. Если новые данные из-

менили стохастическую

еМВ)

'(В)

1

Начало

17) точная оценка параметров модели. Выполняется в соответствии с алгоритмом параметрической оптимизации Марквардта [1], Левенберга-Марк-вардта либо другими методами;

18) модель yt адекватна? Если «да», то переходим к 10-му этапу, если «нет» - к 15-му этапу. Адекватность проверяем в соответствии со значениями статистик, определяемых уравнениями (2), (3);

19) вычисление прогнозов yt(l), l = l,Ly . Осуществляем по уравнению (10). Дисперсию и вероятностные пределы прогнозов выполняем в соответствии с уравнением (19), (25);

20)печать прогнозов yt(l), 1 = 1,1/;

21) ввод фактических значений реализаций процессов yt+i и х^., i = l,N, в t+j момент времени;

22) вычисление остаточных ошибок моделей аЦ-,

4-ь,

i = 1,N или

ф+.(в)

8? (В)

жено на 18-м этапе;

29) окончание алгоритма.

структуру моделей at, то это будет обнару-

Идентификация параметров модели дМ =F/(0)a,(i), i = lN.

Грубая оценка параметров модели

Точная оценка параметров модели

' /

Идентификация периодов и

порядка взятия разностей

SuS2,...,S„s-,

DX,D2,...,D„S.

Преобразование рядов хР, i = lJV, у,

^vJv^.vN« i=m

]- = 1,Ьх ,1 = 1,Киа(,,]- = Ш

выполняется по уравнениям (23), (24);

23) проверка условия: если ■/< 17 , то переходим к 24-му этапу, в противном случае - к 26-му этапу;

24) коррекция прогнозов Уы(1)> 1 = 1 + 1,Ьу. Осуществляется в соответствии с (22);

25)печать скорректированные прогнозов у1+;(1), 1 = 1 + 1,17; '

26) проверка условия «Продолжать ли вычисления?». Если «да», то переходим к 27-му этапу, если «нет» - к 29-му этапу;

27) формирование новых обучающих выборок у1 и х[1), 1 = и соответствующих массивов остаточных

Вычисление прогнозов х^ (/), i = lJV, / = hl.

Вычисление прогнозов j>,(/), / = й; 7 = 1-

Рис. 1. Блок-схема адаптивного алгоритма прогнозирования ВР моделями SARIMAX, расширенными для случая прогнозирования ВР с несколькими сезонными компонентами и с учётом нескольких экзогенных переменных

Рис. 2. Блок-схема адаптивного алгоритма прогнозирования ВР моделями SARIMAX, расширенными для случая прогнозирования ВР с несколькими сезонными компонентами и с учётом нескольких экзогенных переменных

(продолжение)

Алгоритм повышает точность прогнозов ВР в сравнении с алгоритмом для построения моделей SARIMA [4], но требует значительно больших затрат машинного времени и дополнительной информации о внешних факторах.

е(3) = еМ-№(3), t = 1,п;

L3 -1

w((3)=:I 1

4. 2. Математическая модель прогнозирования временных рядов с учётом внешних факторов и с локальной нерегулярностью

Необходимость выполнения предположений о стационарности модели АР и обратимости модели СС приводит к тому, что корни соответствующих характеристических полиномов должны быть вещественными и лежать вне или на единичной окружности, а соответствующие параметры моделей должны находится в диапазоне от -1 до +1 [23]. Процессы, удовлетворяющие этим условиям, будем называть регулярными. В докладе рассматривается более широкий класс взаимосвязанных случайных процессов, как с локальной временной нерегулярностью воздействия внешних факторов, так и с локальной нерегулярностью по частоте, амплитуде и фазе для полигармонических трендов. Показано, что построение адекватных математических моделей таких процессов возможно в классе моделей вида (27) с дополнительным полиномом (В) = 1 -(В) от оператора задержки В, коэффициенты которого определяются как коэффициенты

L-1

линейной рекуррентной формулы (В)а( = ^gaat-i

еК») = £К» -1) - wtn"q), ( = 1, п ;

г(Пзеч) = ^^

--еу(В) 'ф^Дв)а

с(п-еч ) _ ч ч ^ ' „эея . К =^^ ТГЗТ а( ;

у(еч-55А (1) = £ (п-е' )(1) + е(п"е' >(1).

(28)

В мультипликативном виде данная модель запишется следующим образом:

п-ея 95еч '* (В)

УГ-55А' = (В)w(t'> + Я-ея'.( )

Ф5еЧ'+-е,.+ (В) (

П-еЯ -»<]"» (В)

УГ^"^)-(В)е['-1)+--—¿г(В)

Ф5еЯ";-еЧ..+ (В) (

(29)

(30)

метода «Гусеница»-$$А [25], применяемого к остаточным ошибкам а( модели (1). Предложенная модель имеет следующий вид:

SARIMAX-SSA

= 1

N + Г-Н

8(')(В) ("ь' фр+ (в) е(-1 (В)

(27)

где е( - остаточные ошибки модели (27), ( - параметр, называемый длиной окна метода «Гусеница»-SSA [2]. При этом корни характеристического полинома е(-1'(В) =0 могут быть как вещественными, так и комплексными и лежать как снаружи, так и внутри единичной окружности [2, 25].

Вернёмся к рассмотрению модели ЛРФ-SSA+ SARIMA [4]. В более сложных случаях лучшие результаты прогнозирования можно получить, используя последовательный метод SSA, в котором после выделения одной из компонент с некоторым (1 анализируется остаток (снова методом SSA) с возможно уже другой длиной окна (2 и т. д. [2]. Тогда гибридная модель на основе последовательного метода SSA и моделей SARIMA будет иметь вид:

£(1) = У( - wt1), ( = 1,п ;

■I1) ( = 1

(В) (В)..^ (В) угЧ-^"' =

е-еЧ'''Ч3еч..-(В) а8еч„'

Ф5еЧ"'+ч...+ (В) ( ,

('-1

где (В) = Х(В), ' = 1,п5еч - полином от

(31)

)=1

введенного выше оператора задержки В, коэффициен-

ты е^, . = 1,(' -1, ' = 1,п которого равны коэффи-

1 (-1

циентам ЛРФ ^ =Хgj wtJ-', ' = 1,п метода «Гусе-.=1

ница»-SSA; w't'), ' = 1,п - компонента, определяемая на этапе диагонального усреднения метода «Гусени-ца»^^, применяемого к ВР е('-1), ' = 1,п5еч, е(0) = у(, ( = 1, п.

Запишем (31) в следующем виде: ^ е5еч"" (В)

seч-SSA »и__ч ■ '_0зеч 111

У

(В)ф-еЧ'"+-еч...+ (В)

(32)

w « = ! еМ;

e(t2) = e(t1)-w(t2), ( = 1,п;

где

(2-1

w (2)=;£ gj2)wtl'• 1

(В)=е((1-1 (В)е'(2-1 (В)..^ (В), (п. = £(ц-1).

п-еч 1 2 п-еч п-еч ~

Назовём такую модель s:q-SSA.

Тогда гибридная модель (27) с использованием последовательного метода SSA примет следующий вид:

п

п

N Ю(')(В)

ySARIMAX-Seч-SSA = ^ с V )

еч.(В)

8(')(В) (-ь' Фр+(В) '(В)

а- (33) где . х'.1 + ^ + ... + Х'Ч; . , . = 1,г; I -

непересекающиеся подмножества множества индексов {^ 2, d} .

Назовём такую модель SARIMAX - S:q-SSA.

4. 3. Многомерный вариант метода «Гусеница»-88А

В этом разделе рассматривается многомерное обобщение метода «Гусеница»-SSA и приводятся алгоритмы двух вариантов прогнозирования (многомерного продолжения) взаимосвязанных нестационарных ВР.

4. 3. 1. Алгоритм анализа временных рядов многомерным вариантом метода «Гусеница»-88А

Пусть требуется получить прогноз N+1 взаимосвязанных нестационарных ВР Y=y1, у2, ..., уп и Х(')=х(') , х2') , ..., Х^ , ' = произвольных длин П^у и п ('), ' = 1, N соответственно.

1. Вложение

Выбираем длину окна ( и строим траекторную матрицу

X =

х';1

к к

^у Х1

(1)

Х(1) 2

х12) X1

х'2)

к IV

хjN) х(^)

х1

(N) 1 =

]=[

= |У Х(1) X'

(1) Х(2)

X

И]

из векторов вложения

х(') = (х(')

1 <j< К и , К и =п й -1-1, ' = 1,N. Здесь Y - траек-торная матрица ряда_У, Х('), ' = - траекторная матрица рядах('), ' = 1,N.

4. Диагональное усреднение

Матрицы Х1, . = 1,г сгруппированного разложения переводятся в систему новых рядов длины п. Для этого они разбиваются следующим образом

Х^ = [^. Х(1) Х12) ... Х^)]. Далее производится диагональное усреднение каждой из матриц Уг и Х[к), к=1, 2,., N } =_1_г, преобразуя их в ВР ^ и х(к), к=1, 2,., N, . = 1,г соответственно. В результате каждая матрица Хг порождает многомерный

(^ х (1) х !2) ... х м)

ленную аддитивную компоненту исходного ряда

V х(1) х(2) ...

ВР

. = 1,г

восстанов-

4. 3. 2. Прогнозирование временных рядов многомерным вариантом метода «Гусеница»-88А

В теории одномерного метода «Гусеница»-SSA вычисление прогноза ВР выполняется непосредственно по линейной рекуррентной формуле, которую порождает траекторное пространство ряда. Траекторная матрица многомерного ВР несимметричная, поэтому далее рассмотрим два варианта метода «Гусеница»-SSA для прогнозирования многомерных ВР.

Путь имеется N+1 ВР Y, х(1), х(2), ..., х(ю одинаковой длины п. Выберем длину окна 1<(<п.

L-продолжение многомерного ВР

^^ хП1++1 хП2+\ ... хП^1) в пространстве столбцов траекторной матрицы имеет вид:

= Уд^(,

(36)

2. Сингулярное разложение

Сформируем матрицу 8 = ХХТ и произведём сингулярное разложение траекторной матрицы X ВР.

Обозначим:

- Х1, Х2, ..., - собственные числа матрицы 5, взятые в порядке убывания ( > Х2 >...> >0);

- и1, и2, ..., и( - ортонормированная система собственных векторов матрицы S, соответствующих этим собственным числам.

Произведём разложение траекторной матрицы

где

1 -V

(37)

.=1

и1, и2, ..., и, - левые сингулярные вектора сингулярного разложения траекторной матрицы ряда;

ху - вектор, состоящий из первых L - 1 компонент вектора х ;

п. - последняя координата вектора и..

X = X1 +... + X,,,

(34)

где

.=1, 2,.,,; , = тах{^¡Ц >0}.

3. Группировка

Разложение (34) в сгруппированном виде может быть записано следующим образом:

X=X, + Xh

-XI,,

(35)

Уд =

уп-(+2 уп-(+3 • • уп

Х(1) Х(1) • хП1)

Х(2) Х(2) • хП2)

- матрица из последних L-1 значений ВР. ^-продолжение многомерного

^ х1

(1)

ВР

хП2+\ ... хП+1) в пространстве строк тра-

екторной матрицы имеет вид:

^+1 = (^даГ wwт)-1 ^^ ,

(38)

где X * ' - вектор, состоящим из всех компонент вектора X, кроме элементов с номерами, кратными К; V V , ..., V,) - правые сингулярные вектора сингулярного разложения траекторной матрицы ряда;

W =

ny ny • n y d

n(1) Я(21) • п(1) d

n12) п22) • п(2) d

x(N)

n2N)

T(N)

U = (U[<N+1),UJ<N+1),^IUJ<N+1));

ZУ = (Уп-К+2 Уn-K+3 Уп )T

Z(i) = (x(i) x(i) x(i) " ^ ~ ^ n-K+2 n-K+3 • •' An

i=1, 2, ..., N;

Z = (Zy Z(1) Z(2) ••• Z(N)) - вектор последних значений всей совокупности рядов.

4. 4. Гибридные математические модели на основе многомерного варианта метода «Гусеница^-SSA и моделей SARIMAX

В данном разделе предложено, используя формулы L- и K-продолжения формировать передаточные функции моделей (1) и (27).

При использовании формулы L-продолжения матрица YÄ формируется следующим образом:

' Уn-L+2 Уп-L+3 • • Уп 1

x(1) А n+Ь,-L+2 Х(1) Лп+Ь,-L+3 • xn?b,

II Х(2) п+Ь2-L+2 Х(2) n+b2-L+3 Х(2) • Xn+b2

X(N) ^n+bN-L+2 X(N) Än+bN-L+3 • X(N) n+bN У

е;-'(в)

ф++.'(в) * ■

yL-MSSA-SARIMA (j) = yL-MSSA (j) + g (j) ;

-в мультипликативной форме: e;-''(B) ,..

У

L-MSSA-SARIMA .

SL-1 (в)ф++..''(в) *

(40)

где ^ (в) =1 -|А1-1В- •••Дадим обозначение этой модели - SARIMA.

Теперь рассмотрим формулу К-продолжения (38). Сделаем следующие обозначения:

^=(I(N+1)(N+1)"WWT)- WUT;

П-1 (в) = 1 - Ч^В - ^K-2B2 - • - ^BK-1;

(41)

^К-1 (В) = *0+1КК-1)В ^0+1).(К-1)-1В2 ^(1+1)<К-1)-(К-2)ВК 1 ,

1 = Ш.

Тогда выражение для гибридной модели на основе метода К-MSSA и модели SARIMA будет иметь вид: - в аддитивной форме:

gt = У*- У*

K-MSSA * = 1 n ;

_ е\...'"(в)

У K-MSSA-SARIMA (j) = У K-MSSA (j)+ gK-MSSA (j) ;

- в мультипликативной форме:

(42)

K-MSSA-SARIMAX _ i=1

Z^k-1 (B)

n-1 (B)

ф+

(B)

V (B)

a*IV. (43)

T

где - прогнозное значение ВР в момент

времени ]>п. При использовании ^же формулы для К-продолжения векторы Z(i), 1 = 1,К приобретают

следующий вид Z(i) = (хП+b (

x(i)

X (i)

) (¡)-К+2 П+О (¡)-К+3 П+О (¡)

1 =

Прогнозные значения для каждого конкретного ВР х.0 , 1 = 1,К определяются одномерным методом «Гусеница»-SSA [25], моделями SARIMA [1], их комбинированными моделями ЛРФ-SSA-SARIMA, Тренд-SSA-SARIMA [4] или декомпозиционными моделями SSA(Z)-SARШA, SSA-SARIMA, SSA-SARIMA+ SARIMA [4].

Гибридная математическая модель на основе выражения 1-продолжения и модели SARIMA запишется следующим образом:

- в аддитивной форме:

g * = У*- У*

A, * = 1,n;

Обозначим эту модель К-MSSA - SARIMAX. Для модели (27) выражения, использующие для формирования передаточной функции формулы 1-продолжения примут вид: - в аддитивной форме:

е. = У. - УL ^ . = й;

g = е;-'(в) .= _ 1

- ■ •' — N ^ ' * *a

ф++'(в)

gLa-1'(B)

e;

УL-MSSA-SARIMA-SSA (j) = yL-MSSA (j) + * '(j); (44)

в мультипликативной форме: е;*."(в)

,L-MSSA-SARIMA-SSA .

Sl-1 (В)ф++..''(B) gLY'(B)

e''.(45)

T

Дадим обозначение этой модели SARШA-SSA.

Выражения, использующие для формирования передаточной функции модели (27) формулы К-продолжения примут вид: - в аддитивной форме:

у ¥=

(в)

П-1 (в)

= (^К-1К 1-1+^к-2К 1-2+- + ^ 1-(К-1)) +

+£(^(,+1Кк-1)Х(1-ь, ^(1+1)-(К-1)-1Х1-'ь1 -1 —

(1+1){К-1)-(К-2)А 1-Ь, -(К-2)

* 1 = К1 - кК-М55А, 1 = 1,п;

1 = К,п , если строим модель (43) или (47).

вХ"(в) '''

—-га1 ;

"" Ф++-'"(в)

1 ^Г-1'"(в)

е1 ;

KK-MSSA-SARIMA-SSA (}) = уК-МББА (}) + * (}) ;

(46)

3. Вычисляется остаточный ВР е1 следующим образом:

- = у1 _уь^шм, 1 = 1,п, если строим модель (39) или (44);

- ^шм = у1 - к К-MSSA, 1 = 1,п, если строим модель (42) или (46);

в мультипликативной форме:

K-MSSA-SARШAX-SSA I_ 1=1

(в)

П-1 (в)

(в) 1

ф+(в) g•Laí' (в) 1

(47)

е^ = у1 - у 1% 1 = 1, п - L +1, если строим модель (40)

или (45);

- = у1 - у¥, 1 = 1,п - К +1, если строим модель (43) или (47).

4. Идентифицируются модели SARIMA для остаточного ВР

Обозначим эту модель К-MSSA - SARIMAX - SSA.

4. 4. 1. Метод на основе совместного использования многомерного метода «Гусеница»-88А и моделей SARIMAX

Ниже приведен алгоритм метода на основе совместного использования многомерного варианта метода «Гусеница»-SSA и моделей SARIMAX, алгоритм идентификации которых приведен в пункте 4. 1. 5.

1. Если строим модель (39), (42), (44) или (46), то переходим к этапу 3.

Применяя многомерный метод «Гусеница»-SSA к исходным ВР, определяются:

(1) (1)

- полиномы от оператора задержки ю((в), 8;/(в), 1 = 1,4 если строим модель (1) или (27);

- вектор SL и полином 8ь-1 (в), если строим модель (40) или (45);

- вектор ¥ и полиномы ^К-1 (в) и ^К- (в), 1 = 1,4 если строим модель (43) или (47).

2. Определяется ВР у 1 следующим образом:

_ N т(р (в) (1) N ) ) )

- У1 =Х^(Г)(в)ХС =Х(51)у1-1 + 52)у1 -2 + — + 8(,)у1-г, ) +

+х цч-ь,-««„,-1 -—-»(м-

1 >Х(ь, + с,) +1, если строим модель (1) или (27);

1=1

- у[SL ] Ь-1у1-1 + IX ] L-2 У 1-2 + — + [А ] lK(-(L_1) ,

1 = Ь,п , если строим модель (40) или (45);

в*, (в)

е1 = *+ , . а1, если используется модель (1) или (27);

р+ v )

в*,..''(в)

- мssA =—*—, если строим модель (39);

р+ '' ^ '

в*,..''(в)

- *К мssA = _ч—__а1", если строим модель (42);

ФР+'"(в)

- еГь =

s в*,.''(в)

ь - —*-а1'', если строим модель (40);

1 ~Ф++''"(в) в;.'"(в) ,,

1 " ф++'"(в)

в**.. ''(в)

pL-MSSA = д* У / ''

е1 =ф++''"(в)а1 .

а1'', если строим модель (43);

5. Идентифицируем ЛРФ

1

- а1 = *а е1, если строим модель (27) или

-I

а"4', если строим модель (33);

1 с^Ч

*

'(в)

1

1 g*L_l'(B) 1

е1', если строим модель (44) или

1 с*»*

*

'(в)

а"4", если для остаточных ошибок а1'

1

используем последовательный метод «Гусеница^-SSA;

a =-

1

' g*L-l"(B) 1

et'', если строим модель (45) или

a*eq''', если для остаточных ошибок at''

t g^ '''(в)

используем последовательный метод «Гусеница^-SSA;

a =-

1

t g La-l "'(В) 1

et''', если строим модель (46) или

'(В)

aSeqIV , если для остаточных ошибок

at''' используем последовательный метод «Гусени-ца^-SSA;

1

aIV =

t "gLa-v (в) 1

etIV, если строим модель (47) или

t »seq V

(В)

aSeqV , если для остаточных ошибок atIV

используем последовательный метод «Гусеница^-SSA.

6. Вычисляются прогнозы с заданным упреждением l:

- y[-MSSA (l), eL-MSSA (l), l = 1,L , если применяется модель (39);

- yf-MSSA (l), ef-MSSA (l), l = Ц, если применяется модель (42);

- переходим к этапу 6, если используем модели (40) или (43).

7. Синтезируется общий прогноз как сумма частных прогнозов:

- yL-mssa-sarima (l) = yL-MSSA (l) + £L-MSSA (l), если используется модель (39);

- y к-mssa-sarima (l) = yK-MSSA (l) + giS-MSSA (l), если используется модель (42).

Или вычисляется прогноз по модели (40) или (43), предварительно перейдя от операторного к разностному уравнению.

4. 5. Гибридные математические модели и методы декомпозиционного подхода к прогнозированию временных рядов с учётом внешних факторов на основе многомерного варианта метода «Гусеница^-SSA и моделей SARIMAX

По аналогии с алгоритмом декомпозиционного подхода к прогнозированию, приведенного в [4] в работе предложен алгоритм декомпозиционного подхода к прогнозированию ВР на основе совместного использования многомерного метода «Гусеница^-SSA и моделей SARIMAX. Повысить эффективность прогнозирования таких декомпозиционных моделей в большинстве случаев можно, используя вместо моделей SARIMAX искусственные нейронные сети (ИНС) для

прогнозирования каждой из компонент разложения. Схема декомпозиционного подхода к прогнозированию нестационарных ВР с учётом внешних факторов при этом примет следующий вид (рис. 3).

Как показано в [4] в _зависимости от способа группировки матриц X; , j = 1,г можно рассмотреть следующие варианты гибридных моделей декомпозиционного подхода к прогнозированию ВР на основе многомерного метода «Гусеница»-$$А и моделей SARIMAX.

1. На этапе 3 «группировка» многомерного варианта метода «Гусеница»-SSA возьмём параметр г^, т. е. всё множество индексов {1, 2, ..., L} разделим на L подмножеств по одному элементу в каждом следующим образом ; = {j} , j = 1, L . Тогда на этапе 4 «диагональное усреднение» многомерного варианта метода «Гусени-ца»-SSA исходные прогнозируемый и экзогенные ВР будут разложены на следующие суммы:

yt = 1уt° , х« = £хf , j = 1,N.

(48)

Далее декомпозиционный подход к прогнозированию нестационарных ВР на основе совместного использования многомерного варианта метода «Гу-сеница»-SSA и моделей SARIMAX предполагает идентифицировать и строить модели SARIMAX для каждого ВР разложения у, х^ , 1 = , j = 1^.

Для удобства сделаем следующие обозначения:

w(tl) = y(tl), i = 1,L , w\

(L+(j-1)L+i) = x j

= xtj,l), i = 1,L , j = 1,N.

Таким образом, гибридная модель декомпозиционного подхода к прогнозированию ВР на основе многомерного варианта метода «Гусеница»-SSA и моделей SARIMAX можно записать в следующем виде:

ю'

,(j) =

= I г

leM,l*j О

(в)

(в) t-

e?>-*(B) w»

—-a.

Ф

(в)

t

M = {1,2, ..., (N + 1)L };

(N+1)L

yt(l)= I wt)(l).

(49)

2. Теперь рассмотрим данную модель в контексте концепции стандартизованного моделирования нестационарных ВР, приведенного в [4, 26]. На этапе 3 «группировка» многомерного метода «Гусеница»-SSA всё множество индексов {1, 2, ..., L} делится на г непересекающихся подмножеств 11, 12, ..., 1г . Тогда на этапе диагонального усреднения метода «Гусеница»-SSA исходный ВР будет разложен на следующую сумму:

yt = Iytl) , хtj) = Iх j , j = 1,N

(я

(50)

если в непересекающиеся подмножества 11, 12, ..., 1г вошли все элементы из множества {1, 2, ..., L}. Или

I yf , xtj)«I xtj,l), j = 1,N .

(51)

IV

в противном случае. Причём ВР у^, х^'1', 1 = 1,г , } = 1,N могут интерпретироваться как трендовые, квазипериодические или_ шумовые составляющие исходных ВР у и х1, 1 = 1,N соответственно.

( ' (1) /л ^ —

Также

вводим обозначения w^, = у;', 1 = 1,L, »МнМ = хао , 1 = 1,г, 1 = 1,N.

Тогда гибридная модель прогнозирования нестационарных ВР, основанная на концепции стандартизованного моделирования и на многомерном методе «Гусеница»-$$А и моделях SARIMAX для случая (50) примет вид:

(в) е»™* '(В)

^Ч^дв) :»ф »'(в)at

w г га р

-at

M = {1,2, ..., (N + l)r };

(N+l)r

yt (1)=I: t') '(1).

А для случая (51) -

„(»'' =

wW :(J)

''(B)

'' = У

е»; ''(в)

^sw!:):,,.''(B) ф^)+..+''(в)t

W w(') P

M = {l,2,..., (N + l)r };

(N+i)r ю:ос:(1)-'''(B) (,) '' eq.yv (в) v

yt = y-:-wW +_3-aV

1=1 »w■■■'"(в) t-

фр+vy (в)t •

1. Применяя многомерный вариант метода «Гусе-ница»-SSA к исходным прогнозируемому и экзогенным ВР определяются ВР

- у(1), х, 1 = 1,L, 1 = 1^, если строится модель (49);

- у(1), х, 1 = 1,г , 1 = 1,N, если строится модель (52) или (53)

и вводятся обозначения

- »(1) = уI1), 1=Ц, = х^, 1=Ц, 1=Щ

- »(1) = у(1), 1 = й, »[г+°-1)г+1) = х^, 1 = й, 1 = Щ

2. Для каждого ВР

- »(1), 1 = 1,(N +1^, если используется модель (49);

(52)

- »(1), 1 = 1,^ + l)L , если используется модель (52) или (53)

производится структурная и параметрическая идентификация моделей SARIMAX (используя алгоритм пункта 4. 1. 5)

„(» -

<:» (в)

= 1

^-:(')

е^дв)

(53)

Дадим модели (49) обозначение MSSA(I) -SARIMAX, модели (52) - MSSA - SARIMAX, а модели (53) - MSSA - SARIMAX+SARIMA.

Прогнозные значения компонент в данной работе предлагается получать при помощи моделей (1). Однако, для тех же целей можно также использовать одномерный или многомерный вариант метода «Гусеница»-SSA [2, 25], модели SARIMA [1] или SARIMAX, комбинированные модели ЛРФ-SSA-SARIMA, Тренд-SSA-SARIMA [4] или декомпозиционные модели SSA(Z)-SARШA, SSA-SARIMA, SSA-SARIMA+SARIMA [4], а также, синтезированные в данной работе модели I-MSSA - SARIMA, K-MSSA - SARIMAX, I-MSSA - SARIMA - SSA, K-MSSA - SARIMAX - SSA или искусственные нейронные сети.

^С.дв)" t-b:i') ф-j)++(в) t '

: См р: +

M = {1,2,..., (N +1) L} для модели (49);

ю:":<Л (в) е:';л '(в) () iSM'i^l SW|i) :(1) ' (в) :W ф::Ш+ ^

: wM P

M = {1,2,..., (N +1) г} для модели (52);

»:«:<»/'Н e:(j * ''(в)

„(1)'' = У_^_:(i) " » 1 ;

:(l)

SW|» :М.. " (в)

\ к» +■■

t-b m <Г>: +

ф»+ ''(в) t

М = {1,2, ..., (N +1) г} для модели (53).

3. Вычисляются прогнозы с заданным упреждением 1

- 1»(1)(1), 1 = 1,^ +1^, для модели (49);

- (1)(1), 1 = 1,(N +1)г, для модели (52) или (53).

4. Синтезируется общий прогноз как сумма частных прогнозов

4. 5. 1. Метод декомпозиционного подхода к прогнозированию нестационарных временных рядов с учётом внешних факторов на основе совместного использования многомерного варианта метода «Гусе-ница^-SSA и моделей SARIMAX

Приведём алгоритм и схему (рис. 3) метода на основе совместного использования многомерного метода «Гусеница^-SSA и моделей SARIMAX декомпозиционного подхода к прогнозированию нестационарных ВР с учётом экзогенных факторов.

!N+1)L ~ м

yt (1)= У :t (1) для модели (49);

i=1 (N+1)r

yt(1)= У :t''(1) для модели (52) i=1

- вычисляется прогноз по модели (53), предварительно приведя её к разностному выражению.

:

ш

Исходный ВР

1:1

J п —г

Группировка №1 компонент разложения ВР у^

Группировка №г компонент разложения ВР у^

Группировка №1 компонент разложения ВР

Группировка №2 компонент разложения ВР

Группировка №г компонент •и1

разложения ВР

БДШМД-модель или ИНС № 1

Группировка №2 компонент разложения ВР БДШМД-модель или

ИНС № 2

БДШМД-модель или ИНС № £ или № г или № (г+1) в зависимости от способа группировки

БДШМД-модель или

ИНС № (¿+1) или № (г+1) или № (г+2) в

зависимости от способа группировки

БДШМД-модель или

ИНС № (¿+1) или № (г+1) или № (г+2) в

зависимости от способа группировки

БДШМД-модель или

ИНС № 2x1 или № 2хг или № 2х(г+1)

в зависимости от способа группировки

Исходный ВР

N

О

Я г

Разложение ВР

методом «Гусеница»^ А

Группировка №1 компонент разложения ВР

Группировка №2 компонент разложения ВР

БДШМД-модель или ИНС № (Их£+1) или № (Ихг+1)

или № (Мх(г+1)+1) в зависимости от способа группировки

БДШМД-модель или ИНС № (Мх£+2) или № (Мхг+2)

или № (Мх(г+1)+2) в зависимости от способа группировки

БДШМД-модель или ИНС № (И+1)х£ или № (И+1)хг

или № (М+1)х(г+1) в зависимости от способа группировки

Рис. 3. Схема декомпозиционного подхода к прогнозированию нестационарных ВР с учётом внешних факторов

5. Результаты исследований эффективности предложенных моделей прогнозирования временных рядов с учётом внешних факторов

Исследование предлагаемых гибридных моделей на основе многомерного варианта метода «Гусени-ца»-SSA и моделей SARIMAX осуществим сопоставляя их результаты прогнозов с результатами прогнозов, получаемых классическими вероятностными моделями SARIMAХ, обобщёнными на случай прогнозирования ВР с несколькими сезонными составляющими и с учётом нескольких внешних факторов. Реализация рассмотренных моделей производилась в математическом пакете MATLAB R2014a.

Тестирование будем проводить на ВР часовых значений потребления электроэнергии объёма 1008 значений, что соответствует 6 неделям (рис. 4), с учётом изменения температуры воздуха (рис. 5). Обучение моделей будет производиться на выборках данных за 5 недель (840 значений), а тестирование - на данных последней недели. После проверки адекватности модели, прогноз будет выполняться одношаговый и производиться скольжение окна до последнего значения ВР.

Как отмечается в [4], структурная идентификация гибридных математических моделей на основе методов «Гусеница»-SSA и Бокса-Дженкинса является достаточно сложной задачей. Анализ собственных чисел траектор-ной матрицы ВР, автокорреляционной и частной автокорреляционной функций ВР не даёт точных значений структурных параметров наиболее эффективной модели. Поэтому под структурной идентификацией приведенных гибридных моделей будем подразумевать перебор различных вариантов разложения исходных ВР и группировки ВР разложения методом MSSA, идентификацию компонент разложения или остаточной компоненты моделями

SARIMAX, по приведенному в пункте 4. 5. 1 алгоритму и выбор наилучшей по точности прогнозирования синтезированной модели из множества альтернатив. Разложение исходного ВР получается тем более детальным, чем больше длина окна. Наиболее детальное разложение достигается при выборе длины окна, приблизительно равной половине длины ряда и кратной периоду периодической компоненты, за исключением рядов конечного ранга [2]. При одновременном выделении разных компонент ряда (трендовой, периодических, шумовых компонент) все аспекты выбора длины окна должны приниматься во внимание. Таким образом, имеет смысл начинать всегда с выбора длины окна близкой к половине ряда и кратной периоду периодической компоненты, если такая присутствует в ряде и её период известен.

Сравнительный анализ эффективности прогнозирования рассмотренными моделями будем осуществлять при помощи статистики RMSE (Root Mean Squared Error):

1

RMSE = J-£(yt - У

где п - количество вычисленных прогнозов, у( - фактические значения ВР, у( - прогнозные значения ВР.

Тестирование будет проводиться следующим образом. Длина окна метода «Гусеница»-SSA будет фиксированной, но модели одного класса будут отличаться количеством первых компонент разложения, используемых в них. Поэтому на приведенных ниже столбчатых диаграммах на оси абсцисс будем откладывать количество первых компонент разложения, используемых в конкретной модели, а на оси ординат - ошибку прогнозирования, выраженную в статистике RMSE.

Рис. 4. График ВР часовых значений потребления электроэнергии за 6 недель (в МВт)

Рис. 5. График ВР часовых значений изменения температуры воздуха за 6 недель (в С°)

i t=i

..■¡¡Г

Выберем длину окна L=396. В модели будем отбирать количества первых компонент разложения кратные периоду 24. Приведенная ниже столбчатая диаграмма (рис. 6) говорит о снижении ошибки прогнозирования при добавлении к прогнозам, полученным при помощи L-продолжения многомерного варианта метода «Гусеница»-$$А, прогнозов, полученных моделями SARIMAX.

Анализ диаграммы (рис. 7) говорит о снижении ошибки прогнозирования при использовании модели MSSA - SARIMAX+SARIMA в сравнении с прогнозами, полученными моделями MSSA - SARIMAX.

Наименьшую ошибку прогнозирования показала модель K-MSSA - SARIMAX при использовании в модели 312 первых собственных троек SVD-разло-жения траекторной матрицы исходных ВР. Изменяя количество отбираемых в модель первых собственных троек вблизи значения г=312, удаётся найти более точную модель с г=298 (RMSE=8.712). На следующей диаграмме (рис. 9) изображены ошибки прогнозирования RMSE для наилучших моделей из рассмотренных классов моделей. Наименьшую ошибку прогнозирования показала модель K-MSSA - SARIMAX.

В любом случае, всегда стоит попробовать несколько вариантов выбора длины окна. Для рядов со слишком сложной структурой слишком большая длина окна L может привести к нежелательному разложению компонент ряда, в частности, к смешиванию их с другими компонентами ряда. Даже маленькое изменение значения L может уменьшить смешивание компонент разложения и привести к лучшему разделению компонент ряда, т. е. обеспечить переход от слабой к сильной разделимости [2].

На рис. 10, 11 показано, что более точную модель K-MSSA - SARIMAX можно получить, используя длину окна L=72 и количество первых собственных троек SVD-разложения траекторной матрицы исходных ВР метода MSSA равное 46 (RMSE=7.1304). Незначительно по ошибке прогнозирования уступает более простая модель, использующая длину окна L=48 и количество первых собственных троек равное 34. Ошибка прогнозирования такой модели составила RMSE=7.32. Приведём прогнозы наилучшей по эффективност прогнозирования из класса K-MSSA - SARIMAX модели (рис. 12).

Рис. 6. Диаграммы зависимости ошибки прогнозирования RMSE моделей ¿-MSSA, ¿-MSSA — SARIMAX и К'-МЗЗА — SARIMAX от количества используемых в них первых компонент разложения метода MSSA

Рис. 7. Диаграммы зависимости ошибки прогнозирования RMSE моделей MSSA — SARIMAX и MSSA — SARIMAX+SARIMA от количества используемых в них первых компонент разложения метода MSSA

■ -ММА-ЗАШМ А*

295 296 297 290 299 300 301302 303 304 305 306 307 303 309 310 311312 313 314 315 316 31731В 319 320

Рис. 8. Диаграммы зависимости ошибки прогнозирования RMSE моделей Л_-MSSA — SARIMAX от количества используемых в них первых компонент разложения метода MSSA

Рис. 9. Диаграмма минимальных ошибок прогнозирования RMSE для рассмотренных классов моделей (L=396) 100

Рис. 10. Диаграммы зависимости ошибки прогнозирования RMSE моделей — SARIMAX от количества

используемых в них первых компонент разложения методом MSSA для длин окна 24, 48, 72, 96 и 120 соответственно

Рис. 11. Диаграмма минимальных ошибок прогнозирования RMSE для моделей — SARIMAX,

использующих различные длины окон

Рис. 12. График временного ряда потребления электроэнергии и прогнозы, полученные лучшей по эффективности

прогнозирования из рассмотренных моделей

Рассмотрим пример прогнозирования ВР потребления природного газа. График данного ВР среднесуточного объёма газопотребления представлен на рис. 14. График ВР изменения температуры воздуха представлен на рис. 13. Обучение моделей будем проводить на данных за два года (объём 730 значений). После проверки адекватности модели, прогноз будет производиться одношаговый, также как и для предыдущего примера.

Для ВР потребления природного газа, который имеет полиномиальный, полигармонический и стоха-

стический тренд со сложной корреляционной структурой наилучший результат выдала модель из класса MSSA - SARIMAX+SARIMA (рис. 15). Ниже приведём прогнозы наилучшей по эффективности прогнозирования модели из класса MSSA - SARIMAX+SARIMA (рис. 16).

Таким образом, для ВР различной структуры, наиболее эффективный прогноз дают модели из различных классов.

Рис. 13. График ВР среднесуточных значений изменения температуры воздуха (в С°)

Рис. 14. График ВР среднесуточных данных потребления природного газа (в Мм3/ )

сут

Рис. 15. Диаграмма минимальных ошибок прогнозирования RMSE для рассмотренных классов моделей ^=358)

12000 10000 aooo 6000 4000

гооо о

RneMeuunu пал

»^PIII-iiiium

T IfJUI Ul-13

щ IS Г^ со Ol О

Рис. 16. График потребления природного газа и прогнозы, лучшей из рассмотренных моделей

6. Выводы

В работе предложены модели различной сложности для прогнозирования нестационарных ВР с учётом экзогенных факторов, а также методы идентификации этих моделей на основе совместного использования многомерного варианта метода «Гусе-ница^-SSA и моделей SARIMAX, обобщённых для случая прогнозирования ВР с несколькими сезонными компонентами и с учётом нескольких экзогенных переменных. Один класс гибридных математических методов, основанных на многомерном варианте метода «Гусеница^-SSA и моделях SARIMAX, заключается в структурной идентификации передаточной функции модели Бокса-Дженкинса, обобщённой на случай прогнозирования с учётом нескольких экзогенных факторов, используя формулы L- или K-продолжения многомерного варианта метода «Гусеница^-SSA. Декомпозиционный подход к прогнозированию ВР на основе многомерного варианта метода «Гусеница^-SSA и моделей SARIMAX заключается в разложении методом MSSA исходных прогнозируемого и экзогенных ВР на множество ВР с более простой структурой; прогнозировании данных компонент разложения моделями SARIMAX и вычислении общего прогноза, объединяя прогнозы построенных упрощённых моделей. Выбор

полученные наи-

различных параметров многомерного варианта метода «Гусеница^-SSA приводит к созданию, в сочетании с моделями SARIMAX, различных вариантов методов прогнозирования и к синтезу на их основе гибридных математических моделей с различными структурами.

Предложенные модели были протестированы на ВР потребления электроэнергии и природного газа, а их результаты прогнозирования были сравнены с результатами, полученными вероятностными моделями SARIMAX. Экспериментальные результаты показывают высокую эффективность предложенных моделей прогнозирования при выборе подходящих структурных параметров.

Таким образом, в статье разработана методология моделирования и прогнозирования нестационарных ВР, использование которой позволило:

1. Учесть влияние на рассматриваемые ВР произвольного количества основных групп факторов (например, метеорологических, хронологических и организационных).

2. Получить экономичные по структуре и количеству оцениваемых параметров модели, описывающие широкий класс нестационарных ВР с детерминированными (полигармоническими и полиномиальными) или стохастическими трендами и локальной нерегулярностью.

3. Разработать эффективные методы и алгоритмы структурной идентификации моделей, проверки степени их адекватности реальным процессам, вычисления и коррекции прогнозов.

В дальнейшем целесообразно обобщит гибридные математические модели и методы, основанные на многомерном варианте метода «Гусеница^-SSA и моделях VARMAX для прогнозирования взаимосвязанных нестационарных ВР.

Литература

1. Евдокимов, А. Г. Оперативное управление потокораспределением в инженерных сетях [Текст] / А. Г. Евдокимов, А. Д. Тевя-шев. - Х.: Вища школа, 1980. - 144 с.

2. Голяндина, Н. Э. Метод «Гусеница^-SSA: анализ временных рядов [Текст]: учеб. пос. / Н. Э. Голяндина. - СПб.: С.-Петербургский государственный университет, 2004. - 74 с.

3. Варианты метода «Гусеница^-SSA для прогноза многомерных временных рядов [Текст] : труды IV междунар. конф. / Москва "Идентификация систем и задачи управления". - Москва, 2005. - С. 1831-1848.

4. Щелкалин, В. Н. Гибридные модели и методы прогнозирования временных рядов на основе методов «Гусеница^-SSA и Бок-са-Дженкинса [Текст] / В. Н. Щелкалин // Восточно-Европейский журнал передовых технологий. - 2014. - Т. 5, № 4 (71). -С. 43-62. - Режим доступа: http://journals.uran.ua/eejet/article/view/28172/25654 doi: 10.15587/1729-4061.2014.28172

5. Lawrance, A. J. Stochasting modelling of riverflow time series [Text] / A. J. Lawrance, N. T. Kottegoda // Journal of the Royal Statistical Society. Series A (General). - 1977. - Vol. 140, Issue 1. - P. 1-47. doi: 10.2307/2344516

6. 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, 1994. - P. 310-313.

7. Yurekli, K. Application of linear stochastic models to monthly flow data of Kelkit Stream [Text] / K. Yurekli, A. Kurunca, F. Ozturkb // Ecological Modelling. - 2005. - Vol. 183, Issue 1. - P. 67-75. doi: 10.1016/j.ecolmodel.2004.08.001

8. Fraedrich, K. Estimating the dimension of weather and climate attractor [Text] / K. Fraedrich // Journal of the Atmospheric Sciences. - 1986. - Vol. 43. - P. 419-432.

9. Vautard, R. Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series [Text] / R. Vautard, M. Ghil // Physica D: Nonlinear Phenomena. - 1989. - Vol. 35, Issue 3. - P. 395-424. doi: 10.1016/0167-2789(89)90077-8

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

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 // Surveys in Geophysics. - 1996. - Vol. 17, Issue 6. - P. 619-663. doi: 10.1007/bf01931784

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. - P. 320.

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. - Vol. 31. - P. 11721179.

Hassani, H. Forecasting European industrial production with singular spectrum analysis [Text] / H. Hassani, S. Heravi, A. Zhigljavscky // International Journal of Forecasting. - 2009. - Vol. 25, Issue 1. - P.103-118. doi: 10.1016/j.ijforecast.2008.09.007 Щелкалин, В. Н. Трендовый подход прогнозирования временных рядов на основе метода «Гусеница^-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 Power Tech. 2007. - P. 1090-1093. doi: 10.1109/pct.2007.4538467

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 Resources Management. - 2011. - Vol. 25, Issue 11. - P. 2683-2703. doi: 10.1007/s11269-011-9833-y

Щелкалин, В. Н. Декомпозиционный подход прогнозирования временных рядов на основе метода «Гусеница^-SSA [Текст] : Материалы 14-й Межд. науч.-тех. конф. SAIT / В.Н. Щелкалин // УНК "ИПСА" НТУУ "КПИ". - К.: УНК "ИПСА" НТУУ "КПИ", 2012. - С. 260-261.

Щелкалин, В. Н. От идей методов «Гусеница^-SSA и Бокса-Дженкинса до декомпозиционного метода прогнозирования и декомпозиционной ИНС [Текст] / В. Н. Щелкалин // Восточно-Европейский журнал передовых технологий. - 2011. -Т. 4, № 4 (52). - С. 59-69. - Режим доступа: http://www.nbuv.gov.ua/oldJrn/natural/Vejpt/2011_4_4/2011_4_4/59-69.pdf de Souza, R. M. SSA, Hierarchical Clustering and ARIMA Models in Forecasting Reservoir Inflows [Text] / R. M. de Souza, M. Barros, M. Menezes, J. F. M. Pessanha, L. A. T. Junior, K. M. Cassiano, R. C. Souza // The 32nd Annual International Symposium on Forecasting ISF - 2012, 2012. - P. 11-25.

de Souza, R. M. Streamflow Forecasting Using Singular Spectrum Analysis, ARIMA Model And Artificial Neural Network [Text] / R. M. de Souza, M. Menezes, R. C. Souza, J. F. M. Pessanha, L. A. T. Junior, K. M. Cassiano // The 32nd Annual International Symposium on Forecasting ISF - 2012, Boston, 2012. - P. 33-47.

Cassiano, K. ARIMA Model, Neural Networks and SSA in the Short Term Electric Load Forecast [Text] / K. Cassiano, M. Menezes, L. A. Junior, J. Pessanha, Raf. Souza, R. Souza // Economic Forecasting - Past, Present and Future: International Symposium of Forecasting ISF - 2014, Rotterdam, 2014. - P. 1 - 23.

Бокс, Дж. Анализ временных рядов. Прогноз и управление. Вып. I. Пер. с англ. [Текст] / Дж. Бокс, Г. Дженкинс. - М.: Мир, 1974. - 406 с.

Евдокимов, А. Г. Об одном алгоритме выявления скрытой периодичности [Текст] / А. Г. Евдокимов, А. Д. Тевяшев // Автоматизированные системы управления и приборы автоматики. - 1977. - Вып. 43. - С. 60-64.

Голяндина, Н. Э. Метод «Гусеница^-SSA: прогноз временных рядов [Текст]: учеб. пособие / Н. Э. Голяндина - СПб. : С.-Петербургский государственный университет, 2004. - 52 с.

Седов, А. В. Моделирование объектов с дискретно-распределёнными параметрами: декомпозиционный подход [Текст] / А. В. Седов. - Южный научный центр РАН. - М. : Наука, 2010. - 438 с.

Чумак, О. В. Энтропии и фракталы в анализе данных [Текст] / О. В. Чумак. - НИЦ «Регулярная и хаотическая динамика». - М. Ижевск: Институт компьютерных исследований, 2011. - 164 с.