Прогнозування часових рядів методом SSA (приклад)
Матеріал із MachineLearning.
SSA (Singular Spectrum Analysis, "Гусениця")- метод аналізу та прогнозу часових рядів. Базовий варіант методу полягає в:
1) перетворенні одновимірного ряду на багатовимірний за допомогою однопараметричної зсувної процедури (звідси і назва "Гусениця");
2) дослідженні отриманої багатовимірної траєкторії за допомогою аналізу основних компонентів (сингулярного розкладання);
3) відновлення (апроксимації) ряду за вибраними головними компонентами.
Таким чином, результатом застосування методу є розкладання часового ряду на прості компоненти: повільні тренди, сезонні та інші періодичні або коливальні складові, а також компоненти шумові. Отримане розкладання може бути основою прогнозування як низки, і його окремих складових. "Гусениця" допускає природне узагальнення на багатовимірні часові ряди, і навіть у разі аналізу зображень. У статті розглянемо варіант алгоритму, призначений для аналізу багатовимірного часового ряду. Алгоритм містить два вхідні параметри довжину гусениці та її компонент, вибір яких істотно впливає результат роботи алгоритму.
Зміст
Постановка задачі
Спостерігається система функцій дискретного аргументу < де k = 1, . s>. s – число тимчасових рядів, k – номер ряду, N – довжина тимчасового ряду, i – номер відліку. Потрібно розкласти ряд у суму компонентів (використовуючи метод головних компонентів, див. опис алгоритму), інтерпретувати кожну компоненту, і побудувати продовження ряду вибраних компонентів.
Опис алгоритму
Побудова матриці спостережень
Розглянемо спочатку одномірний часовий ряд Виберемо n таке, що час життя багатовимірної гусениці. Нехай -довжина гусениці. Побудуємо послідовність з n векторів у наступному виді:
Зватимемо нецентрованою матрицею спостережень, породженою гусеницею з часом життя n.
У разі багатовимірного часового ряду матрицею спостереження називається стовпець із матриць спостережень, що відповідають кожній з компонентів.
Проведений надалі аналіз основних компонентів може проводитися як по центрованій, так і по нецентрованій вибірках. Для спрощення викладок розглянемо найпростіший нецентрований варіант.
Аналіз основних компонентів
Розглянемо підступну матрицю отриманої вибірки:
Виконаємо її svd-розкладання:
де - діагональна матриця власних чисел, - ортогональна матриця власних векторів.
Далі розглянемо систему основних компонент:
Після проведення аналізу основних компонентів зазвичай передбачається проведення операції відновлення вихідної матриці спостережень по деякому піднабору основних компонентів, тобто для обчислюється матриця .
Далі відновлюються вихідні послідовності. В одновимірному випадку i-а компонента відновленого ряду є середнє значення по i-ої діагоналі відновленої матриці спостережень:
У багатовимірному випадку усереднення проводиться з урахуванням того, що матриця спостережень складається з підматриць, що відповідають кожній компоненті ряду:
Числовий ряд називається продовженням ряду, якщо породжена ним при гусеничній обробці вибірка лежить у тій же гіперплощині, що й у вихідного ряду. Нехай у нас є певний набір обраних головних компонентів Визначимо
Тоді прогнозовані значення системи у точці обчислюються за такою формулою:
Чисельний експеримент
Модельні дані
Розглянемо тривимірний часовий ряд, заданийформулами:
Крім тренду (лінійна, квадратична функції та логарифм відповідно) до рядів додано дві періодичні складові (синусоїди з різними амплітудами та фазами). Також до тимчасових рядів додано шум із середньоквадратичним відхиленням 0,5:
Досліджуємо основні компоненти даного часового ряду, використовуючи гусеницю довжини 50.
Логарифми перших двадцяти власних чисел коварійної матриці:
З графіка видно, що змістовний зміст несуть, мабуть, перші сім основних компонент, інші, швидше за все, шум.
Парі своїх власних чисел відповідають основні компоненти, що відповідають одній частоті. Відкладемо одному графіку основні вектори, відповідальні 3-му і 4-му власним числам. Покажемо, що вони корелюють та відповідають періодичної складової з періодом 11:
Відновивши тимчасовий ряд (на малюнку зображений перший із трьох) за цими двома головними компонентами, переконуємося в цьому:
Аналогічно, 5-му та 6-му власним числам відповідає періодична складова з періодом 8.
1-а, 2-а та 7-а головна компоненти відповідають за тренд:
Відновимо ряд за іншими основними компонентами, припускаючи їх шумовими:
Спрогнозуємо тимчасовий ряд на 50 одиниць вперед за першими семи головними компонентами:
Бачимо, що прогноз цілком адекватний.
Порушення періодичності низки
Подивимося, як відреагує алгоритм порушення періодичності часового ряду. У другій половині часового ряду змінимо частоту, амплітуду і фазу однієї з періодиків:
У цьому випадку 2-а та 3-та головна компонента буде відповідати тій періодиці, яку ми не змінили.
5-а та 6-а відповідає першійполовині зміненої компоненти:
7-а та 8-а - другій її половині:
1-а, 2-я та 9-а відповідають тренду:
Отже, алгоритм стійкий до порушень періодичності - різним ділянкам часового ряду відповідають різні основні компоненти.
Подивимося поведінку алгоритму у разі викидів. Ми додали викид у середньому у кожну 50-у точку часового ряду. Перше власне число в цьому випадку на кілька порядків більше за інших, графік логарифмів власних чисел з 2-го по 20-е представлений нижче:
Відновлюючи ряд перших 16-ти головних компонентів, бачимо, що алгоритм працює некоректно:
Реальні дані
Розглянемо щомісячні дані про покриття півкуль Сонця плямами з 1876 року.
Для північної півкулі:
Використовуємо гусеницю довжини 250.
2-а та 3-я компоненти відповідають 12-річним сонячним циклам:
Автор статті не зміг інтерпретувати решту основних компонентів. Проте суб'єктивно перші шість головних цілком адекватно характеризують сонячну активність і може бути використані прогнозу.