4.2 Базовый алгоритм сингулярного спектрального анализа

Построить математическую модель для анализа внутренней структуры и прогнозирования временного ряда можно с помощью алгоритма ССА, основанного на представлении исходного ряда в виде многомерного, с последующим вычислением собственных векторов и собственных значений ковариационной матрицы исходных данных.

ССА зародился в 70-х годах прошлого столетия из теории динамических систем. В зарубежной литературе метод наиболее известен под названием Singular Spectrum Analysis (SSA) [149, 150, 180, 191, 203, 208, 209]. В его основе лежит трансформация временного ряда в матрицу и ее сингулярное разложение [44, 135]. После идентификации компонент сингулярного разложения происходит их группировка, приводящая к разложению исходного ряда на аддитивные компоненты, такие как тренд, периодики и шум [113, 124, 146, 155, 160, 167, 185, 190, 192, 197]. При этом данный метод не требует стационарности ряда, знания заранее определенной модели тренда, а также сведений о наличии в ряде периодических составляющих и их периодах. Также с помощью данного метода можно определить модель тренда и использовать это знание для дальнейшей обработки ряда уже с известной моделью.

Математической основой метода ССА является сингулярное разложение [86, 144]. Для успешного применения метода ССА следует последовательно пройти несколько шагов. Пусть имеется временной ряд, состоящий из действительных положительных элементов, длины , причем . Приведем описание алгоритма ССА [11, 15, 16, 17, 103, 104, 105]. Алгоритм состоит из выполнения четырех этапов: вложение, сингулярное раз­ложение, группировка и диагональное усреднение. Первые два в совокуп­ности называются разложением, последние два – восстановлением. Основным параметром алгоритма служит так называемая длина окна , такая что . Ре­зультатом алгоритма является разбиение временного ряда на аддитивные составляющие.

Первым этапом является процедура вложения, которая переводит исходный временной ряд в по­следовательность многомерных векторов.

Пусть L — некоторое целое число (длина окна), . Процедура вложения образует так называемых векторов вложения

(4.6)

где имеющих размерность .

L – траекторная матрица исходного временного ряда вида:

(4.7)

состоит из векторов вложения (4.6) в качестве столбцов. Другими словами, траекторная матрица – это матрица, которую можно представить как:

.

(4.8)

Очевидно, что и матрица X имеет одинаковые эле­менты на «диагоналях» . Таким образом, траекторная матрица является ганкелевой. Существует взаимно-однозначное со­ответствие между ганкелевыми матрицами размерности и рядами длины .

Из-за того, что нет общих рекомендаций по выбору ширины окна, параметр L зависит от решаемой задачи и предварительной информации, известной о ряде. Например, для выделения тренда рекомендуется выбирать ширину окна не слишком большой. С другой стороны, для выделения гармонических колебаний рекомендуется большая ширина окна.

Далее проводится сингулярное разложение матрицы A. Матрица симметричная и неотрицательно определенная, а значит ее собственные числа вещественны и неотрицательны. Представленные в упорядоченном виде собственные числа называют сингулярными значениями матрицы A. Пусть – соответствующие им ортонормированные собственные вектора, а - порядок сингулярного разложения. Обозначим:

(4.9)

где .

Тогда сингулярным разложением матрицы A называется ее представление в виде суммы элементарных матриц:

.

(4.10)

Каждая из матриц имеет ранг, равный единице. Поэтому их можно назвать элементарными матрицами. Вектор называют k-м левым сингулярным вектором или просто k-м собственным вектором, вектор – правым сингулярным вектором. Набор будем называть k-ой собственной тройкой.

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

Процедура группировки формально одинакова для всех разновидностей ССА. На основе разложения (4.10) процедура группировки делит все множество индексов на непересекающихся подмножеств .

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

.

(4.11)

Процедура выбора множеств и называется группировкой собственных троек.

Этап диагонального усреднения. На последнем шаге базового алгоритма каждая матрица сгруппированного разложения переводится в новый ряд длины . Для произвольной матрицы X процедура приведения ее к ганкелевому виду и последующему преобразованию в ряд G выражается следующим образом. Пусть – матрица размера с элементами , , . Полагая, что , , , если и в остальных случаях, диагональное усреднение переводит матрицу в ряд по формуле:

(4.12)

Это выражение соответствует усреднению элементов матрицы вдоль побочных диагоналей : выбор дает , для получаем и так далее. Применив диагональное усреднение к матрицам, полученным на этапе группировки, приходим к разложению исходного ряда в сумму различных рядов [16, 17].

Сравнение метода ССА целесообразно проводить с двумя другими распространенными методами структурного анализа временных рядов – Фурье- и Вейвлет-анализом [32, 33, 89]. Само сравнение можно провести с помощью эксперимента методом имитационного моделирования [141]. Пусть:

(4.13)

где – экспоненциальный тренд, – периодическая компонента, – случайная составляющая, полученная как нормально распределенный белый шум с нулевым математическим ожиданием и дисперсией .

В качестве задачи сформулируем необходимость определения и выделения экспоненциального тренда и окрашенного шума из исходного временного ряда [85, 109]. График смоделированных значений ряда представлен на рис. 4.1.

Рис. 4.1. Временной ряд

Спектр Фурье временного ряда представлен на рис. 4.2. Вейвлет спектр представлен на рис. 4.3. На рис. 4.4 представлены первые 10 главных компонент ССА – разложения сигнала, с длиной окна 40. Таким образом, наличие экспоненциального тренда и окрашенного шума не позволяет выявить реально присутствующие в ряде периодические составляющие с помощью Фурье-анализа или дискретизированного вейвлет-преобразования, в отличие от метода ССА, справившегося с поставленной задачей.

Рис. 4.2. Фурье-спектр

Рис. 4.3. Вейвлет-спектр

Рис. 4.4. ССА разложение y[n]

Обобщая полученные результаты, произведем сравнение методов структурного анализа временных рядов (таблица 4.1 и 4.2) [110, 116, 196, 206, 207].

Таблица 4.1

Особенности методов структурного анализа временных рядов

Метод	Вид оконного преобразования	Показания к использованию	Противопоказания к использованию
Фурье-преобразование	В порядке убывания уровня боковых лепестков: •Прямоугольное (-13 дБ). •Ханна (-31.5 дБ). •Хемминга (-42 дБ). •Блэкмана (-58 дБ). •Блэкмана-Харриса (-92 дБ). Применяются из-за невозможности одновременного обеспечения хорошего разрешения по времени и по частоте. Чем уже окно, тем выше разрешение по времени и ниже разрешение по частоте.	Быстрый анализ сигналов несложной формы с ярко выраженными гармоническими составляющими, например кардиосигналов; детрендированных слабозашумленных временных рядов.	- Отсутствие временной развертки. - Невозможность интерпретации сигналов с ярко выраженными трендами (требование стационарности и эргодичности). - Необходимость соблюдения фазовых сдвигов. - Трудность восстановления исходного сигнала.

(Продолжение таблицы 4.1)

Вейвлет-преобразование

- МНАТ (хорошо приспособлен для анализа сложных сигналов). - Morlet (высокая угловая избирательность базиса). - Paul (часто применяется в квантовой механике). - LMB (служит основой для построения ортогональных дискретных базисов по процедуре Малла). - Daubechies (биортогональный, имеющий пару, используемую для получения реконструкционной формулы). Временное разрешение увеличивается с частотой (частотное снижается).

- Быстрый анализ сигналов с частотно-временной разверткой. - Анализ сигналов с переменной частотой. - Широко применяется при анализе изотропных полей, сжатии данных.

- Большая по сравнению с Фурье преобразованием сложность вычислений. - Невозможность интерпретации сигналов с ярко выраженными трендами.

(Продолжение таблицы 4.1)

Сингулярный спектральный анализ

Используется только оконное преобразование длины L - «свёртка» временного ряда в матрицу, содержащую фрагменты временного ряда, полученные с некоторым сдвигом, для перехода от одномерного временного ряда в многомерный ряд и последующего применения к полученному многомерному временному ряду метода главных компонент.

- Анализ сигналов сложной формы. - Прогнозирование. - Интерпретация внутренней структуры временного ряда. - Не требует стационарности. Широко применяется при анализе и прогнозировании временных рядов наблюдений – биологических, экономических, медицинских, технических, например кардиограммы, анализ вариабельности сердечного ритма, анализ потока технических инцидентов в системе газораспределения.

- Большая сложность вычислений. - Сложность анализа сигналов переменной частоты. - Плохая разделимость компонент при некоторых значениях длины окна. - Отсутствие универсальных алгоритмов группировки компонент и выбора длины окна.

Таблица 4.2

Характеристика методов структурного анализа временных рядов

Метод	Преимущества	Недостатки
Анализ Фурье	Раскрывается физическая природа сигнала.	1. Требование стационарности процесса. 2. Требование эргодичности процесса. 3. Необходимость соблюдения фазовых сдвигов. 4. Трудность восстановления исходного сигнала (при прогнозировании).
Вейвлет-анализ	1. Возможность физической интерпретации за счет частотно-временной локализации сигнала. 2. Двумерное представление по системам ортогональных функций.	1. Низкая скорость обработки данных. 2. Жестко фиксируемые параметры. 3. Однозначность преобразования.
ССА	1. Способность осуществлять прогноз согласованный с различными вариантами разложения. 2. Компоненты разложения определяются параметром разложения (длиной окна). 3. Набор сингулярных составляющих порождается исследуемой функцией. 4. Возможность управляемого восстановления процесса.	1. Отсутствие универсального алгоритма группировки компонент. 2. Получение искаженных компонент разложения.