SSA (метод)

SSA (Singular spectrum analysis или Анализ сингулярного спектра) — метод анализа временных рядов, основанный на преобразовании одномерного временного ряда в многомерный ряд с последующим применением к полученному многомерному временному ряду метода главных компонент.

Способ преобразования одномерного ряда в многомерный представляет собой «свёртку» временного ряда в матрицу, содержащую фрагменты временного ряда, полученные с некоторым сдвигом. Общий вид сдвиговой процедуры напоминает «гусеницу», поэтому сам метод нередко так и называют — «Гусеница»: длина фрагмента называется длиной «гусеницы», а величина сдвига одного фрагмента относительно другого шагом «гусеницы». Обычно используется шаг 1.

Singular spectrum analysis (SSA) сочетает в себе элементы классического анализа временных рядов, многомерной статистики, многомерной геометрии, динамических систем и обработки сигналов. К источникам происхождения SSA можно отнести Метод главных компонент и классическую теорему Карунена-Лоэва для спектрального разложения временных рядов и цифровых изображений.

Диапазон областей знаний, где SSA может быть применён, очень широк: климатология, океанология, геофизика, техника, обработка изображений, медицина, эконометрика и многие другие. Поэтому в практических приложениях используются различные модификации SSA. Можно выделить два главных направления, это SSA как универсальный метод (Golyandina et al, 2001) для решения задач общего назначения, таких как выделение тренда, обнаружение периодичностей, корректировка на сезонность, сглаживание, подавление шума, а также SSA для спектрального анализа стационарных временных рядов (Vautard and Ghil, 1989), имеющий большое число приложений в тех областях, где такие ряды наблюдаются, в частности, в климатологии.

SSA как метод анализа ряда

SSA может быть использован без предварительного задания модели ряда для анализа произвольных, в том числе, нестационарных, рядов. Основная цель SSA — разложить ряд в сумму интерпретируемых компонент, таких как тренд, периодические компоненты, шум. При этом знание параметрической формы этих компонент не требуется.

Рассмотрим вещественно-значный ряд $\mathbb {X} =(x_{1},\ldots ,x_{N})$ длины $N$. Пусть $L$ $\ (1<L<N)$ — некоторое целое число, называемое длина окна, и $K=N-L+1$.

Базовый алгоритм SSA

Шаг 1: Вложение.

Строится $L\!\times \!K$ траекторная матрица ряда $\mathbb {X}$ следующим образом:

$\mathbf {X} =[X_{1}:\ldots :X_{K}]=(x_{ij})_{i,j=1}^{L,K}={\begin{bmatrix}x_{1}&x_{2}&x_{3}&\ldots &x_{K}\\x_{2}&x_{3}&x_{4}&\ldots &x_{K+1}\\x_{3}&x_{4}&x_{5}&\ldots &x_{K+2}\\\vdots &\vdots &\vdots &\ddots &\vdots \\x_{L}&x_{L+1}&x_{L+2}&\ldots &x_{N}\\\end{bmatrix}},$

где $X_{i}=(x_{i},\ldots ,x_{i+L-1})^{\mathrm {T} }\;\quad (1\leq i\leq K)$ — вектора вложения длины $L$. Матрица $\mathbf {X}$ является ганкелевой, то есть $\mathbf {X}$ имеет одинаковые элементы $x_{ij}$ на анти-диагоналях $i+j=\,{\rm {const}}$.

Шаг 2: Сингулярное разложение (SVD).

Производится сингулярное разложение (SVD) траекторной матрицы $\mathbf {X}$. Положим $\mathbf {S} =\mathbf {X} \mathbf {X} ^{\mathrm {T} }$ и обозначим $\lambda _{1},\ldots ,\lambda _{L}$ собственные числа $\mathbf {S}$, взятые в невозрастающем порядке ($\lambda _{1}\geq \ldots \geq \lambda _{L}\geq 0$), и $U_{1},\ldots ,U_{L}$ ортонормированную систему собственных векторов матрицы $\mathbf {S}$, соответствующих собственным числам.

Положим $d=\mathop {\mathrm {rank} } \mathbf {X} =\max\{i:\lambda _{i}>0\}$ (заметим, что в реальности, как правило, $d=L$) и $V_{i}=\mathbf {X} ^{\mathrm {T} }U_{i}/{\sqrt {\lambda _{i}}}$ $(i=1,\ldots ,d)$. В этих обозначениях сингулярное разложение траекторной матрицы $\mathbf {X}$ может быть записано как

$\mathbf {X} =\mathbf {X} _{1}+\ldots +\mathbf {X} _{d},$

где матрицы $\mathbf {X} _{i}={\sqrt {\lambda _{i}}}U_{i}V_{i}^{\mathrm {T} }$ имеют ранг 1 и называются элементарными матрицами. Набор $({\sqrt {\lambda _{i}}},U_{i},V_{i})$ называется $i$-й собственной тройкой (коротко — ET от eigentriple) сингулярного разложения. Векторы $U_{i}$ и $V_{i}$ называются левыми и правыми, соответственно, сингулярными векторами матрицы $\mathbf {X}$, числа ${\sqrt {\lambda _{i}}}$ — сингулярные числа (они составляют сингулярный спектр $\mathbf {X}$, что и дало название Singular Spectrum Analysis методу), векторы ${\sqrt {\lambda _{i}}}V_{i}=\mathbf {X} ^{\mathrm {T} }U_{i}$, по аналогии с анализом главных векторов, называются векторами главных компонент.

Шаг 3: Группировка собственных троек.

Множество всех индексов $\{1,\ldots ,d\}$ разбивается на $m$ непересекающихся подмножеств $I_{1},\ldots ,I_{m}$.

Пусть $I=\{i_{1},\ldots ,i_{p}\}$. Тогда результирующая матрица $\mathbf {X} _{I}$, соответствующая группе $I$, определяется как $\mathbf {X} _{I}=\mathbf {X} _{i_{1}}+\ldots +\mathbf {X} _{i_{p}}$. Результирующие матрицы вычисляются по группам $I=I_{1},\ldots ,I_{m}$ и сгруппированное SVD разложение матрицы $\mathbf {X}$ может быть записано как

$\mathbf {X} =\mathbf {X} _{I_{1}}+\ldots +\mathbf {X} _{I_{m}}.$

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

Каждая матрица $\mathbf {X} _{I_{j}}$ сгруппированного разложения ганкелизуется (усредняется по анти-диагоналям) и затем полученная ганкелева матрица трансформируется в новый временной ряд длины $N$ на основе взаимно-однозначного соответствия между ганкелевыми матрицами и временными рядами. Диагональное усреднение, применённое к каждой результирующей матрице $\mathbf {X} _{I_{k}}$, производит восстановленные ряды ${\widetilde {\mathbb {X} }}^{(k)}=({\widetilde {x}}_{1}^{(k)},\ldots ,{\widetilde {x}}_{N}^{(k)})$. Таким образом, исходный ряд $x_{1},\ldots ,x_{N}$ раскладывается в сумму $m$ восстановленных рядов:

$x_{n}=\sum \limits _{k=1}^{m}{\widetilde {x}}_{n}^{(k)}\ \ (n=1,2,\ldots ,N).$

Данное разложение является главным результатом алгоритма SSA для анализа временного ряда. Это разложение имеет смысл. если каждая из его компонент может быть интерпретируема как либо тренд, либо колебания (периодики), либо шум.

Метод имеет такие модификации как SSA с однократным центированием и SSA с двойным центрированием. Последний вариант хорошо работает при наличии линейного тренда.

Теория SSA

Теория SSA отвечает на следующие вопросы: (a) какие составлящие временного ряда могут быть разделены SSA и (b) как выбрать длину окна и провести правильную группировку, чтобы выделить нужную компоненту. Основные теоретические результаты содержатся в Golyandina et al (2001, Ch. 1 and 6).

Тренд (определяемый как медленно-меняющаяся компонента ряда), периодические компоненты и шум асимптотически разделимы SSA при $N\rightarrow \infty$. На практике $N$ фиксировано и речь идёт о приближенной разделимости компонент временного ряда. Определить, есть ли приближенная разделимость, можно с помощью ряда индикаторов, см. Golyandina et al (2001, Ch. 1). Длина окна $L$ определяет разрешение метода: большие значения $L$ (но не более половины длины ряда) обеспечивают наиболее подробное разделение на элементарные компоненты и, как следствие, лучшую разделимость. В некотором смысле, длина окна $L$ определяет разрешимость метода, в частности, $L$ соответствует максимальному периоду, который может быть обнаружен с такой длиной окна. Тренд может быть выделен группировкой собственных строек с медленно-меняющимися собственными векторами. Синусоиде с частотой меньше 0.5 соответствует пара синусообразных собственных вектора с той же частотой и разницей в фазах, примерно равной $\pi /2$.

Разделение двух временных рядов может быть формализовано как выделение одного ряда в присутствии возмущения другим рядом. Применение теории возмущений в SSA рассматривается в работе Nekrutkin (2010).

SSA для стационарных рядов

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

Основное отличие алгоритма заключается в использовании лаг-ковариационной матрицы $\mathbf {C}$ вместо $\mathbf {S} =\mathbf {X} \mathbf {X} ^{\mathrm {T} }$. Матрица ${\mathbf {C} }$ может быть оценена непосредственно по исходному ряду как теплицева матрица с константами на диагоналях (Vautard and Ghil, 1989):

$c_{ij}={\frac {1}{N-|i-j|}}\sum _{t=1}^{N-|i-j|}x_{t}x_{t+|i-j|}.$

Использование матрицы $\mathbf {S}$ в базовом методе SSA было предложено в Broomhead and King (1986a,1986b) и поэтому Базовый вариант SSA иногда называют `BK-SSA’. Использование матрицы $\mathbf {C}$ было введено в Vautard and Ghil (1987), что дало методу имя `VG-SSA’ (другое название этой модификации — ‘Теплицев SSA’ (Golyandina et al, 2001, Sect. 1.7)).

Прогнозирование и другие расширения SSA

Прогнозирование. В случае $x_{n}=s_{n}+e_{n}$ на основе SSA разложения можно построить оценку сигнального пространства и получить оценку коэффициентов линейного рекуррентного соотношения, управляющего сигналом, то есть $s_{n}=\sum _{k=1}^{L-1}a_{k}s_{n-k}$. Сигнальное подпространство и полученное линейное рекуррентное соотношение служат основой для алгоритмов SSA прогнозирования, в частности, для рекуррентного и векторного прогнозирования.

Одновременный анализ системы временных рядов (MSSA, многомерный SSA). Если на этапе вложения построить траекторную матрицу системы временных рядов, состыковав вместе траекторные матрицы одномерных рядов, то мы получим метод MSSA, позволяющий строить одновременное разложение сразу нескольких рядов. MSSA лучше применения SSA к рядам по отдельности, если ряды имеют похожую структуру. Можно также рассмотреть одновременный прогноз системы рядов.

Анализ цифровых изображений (2D-SSA). В этом случае аналог траекторной матрицы строится специальным образом с помощью скользящего двумерного окна размера $L_{x}\times L_{y}$.

Заполнение пропусков во временных рядах Существует два метода заполнения пропусков. В алгоритме, описанном в Kondrashov and Ghil (2006), заполнение пропусков проводится на основе итерационной процедуры, когда на каждой итерации на места пропусков ставятся значения, полученный на предыдущей итерации. В Golyandina and Osipov (2007) используется идея восстановления пропущенных компонент векторов из заданного подпространства. Соответственно, рекуррентный и векторный прогноз являются частным случаем заполнения пропусков, если задать пропуски на местах прогнозируемых значений.

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

Соотношение между SSA и другими методами

SSA и Авторегрессия. Типичная модель ряда, рассматриваемая в SSA, это $x_{n}=s_{n}+e_{n}$, где $s_{n}=\sum _{k=1}^{r}b_{k}s_{n-k}$ (сигнал, управляемый линейным рекуррентным соотношением) и $e_{n}$ — шум. Модель авторегрессии (AR) имеет вид < $x_{n}=\sum _{k=1}^{r}b_{k}x_{n-k}+e_{n}$. В первом случае (SSA) шум добавляется ко всему сигналу, а во втором (AR) — на каждом шаге. Несмотря на то, что модели выглядят похожими, SSA рассматривает авторегрессию только как модель для шума.

SSA и Разложение Фурье. В отличие от анализа Фурье, где рассматривается фиксированный базис из синусов и косинусов, SSA использует адаптивный базис, порождаемый самим рядом. В результате, модель ряда, лежащая в основе SSA, более общая, и SSA может выделять амплитудно-модулированные синусы и косинусы с частотами, отличающимися от $k/N$. Соответственно, методы, основанные на подпространстве сигнала, оцениваемом SSA, могут оценивать частоты с более высоким разрешением, чем спектральный Фурье анализ.

SSA и методы, основанные на подпространстве сигнала. SSA можно рассматривать и как метод обработки сигналов, в основе которого лежит оценивание сигнального подпространства, так как оценка сигнального подпространства размерности $r$ может быть получена в рамках SSA как $\mathop {\mathrm {span} } (U_{1},\ldots ,U_{r})$.

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

SSA и линейные фильтры. Восстановление составляющей ряда с помощью SSA может рассматриваться как адаптивная линейная фильтрация. Если длина окна $L$ мала, то каждый собственный вектор $U_{i}=(u_{1},\ldots ,u_{L})^{\mathrm {T} }$ порождает линейный фильтр ширины $2L-1$, который даёт восстановление середины ряда ${\widetilde {x}}_{s}$, $L\leq s\leq K$. Фильтрация не является причинной. Однако, так называемый Causal SSA может быть рассмотрен как аналог причинного фильтра (Golyandina and Zhigljavsky 2013, Sect. 2.9).

Краткая история метода

Первая публикация, которая может рассматриваться как один из источников происхождения SSA и методов обработки сигналов, основанной на оценивании сигнального подпространства, относится к восемнадцатому веку (метод Прони).

Следующая, относящаяся к анализу динамики численности животных, Ефимов, Галактионов (1983). Первая в мире монография на эту тему — Ефимов и др. (1988).

Broomhead and King (1986a, b) и Fraedrich (1986) сформулировали алгоритм SSA для использования его в контексте нелинейной динамики для реконструкции аттракторов.

Ghil, Vautard и их коллеги (Vautard and Ghil, 1989; Ghil and Vautard, 1991; Vautard et al., 1992) заметили аналогию между траекторной матрицей из работ Broomhead и King, с одной стороны, и разложением Карунена-Лоэва (Метод главных компонент во временной области), с другой. Таким образом, SSA стал использоваться как метод для анализа временных рядов независимо от реконструкции аттракторов, включая и те случаи, когда последнее не имеет смысла.

Можно отдельно отметить методологию так называемого метода ‘Гусеница’, см. Данилов и Жиглявский (Ред.) (1997) и Golyandina et al (2001). Эта методология является версией SSA, которая исходно была разработана в СССР независимо от зарубежных работ. Основным отличием методологии «Гусеницы»-SSA является то, что метод разрабатывается для анализа рядов общего вида, основной акцент ставится на изучение теоретических свойств метода. Основным теоретическим понятием является понятие разделимости рядов. В частности, требование разделимости рядов накладывает свои ограничения на выбор параметров.

В настоящее время есть несколько десятков статей с методологическими аспектами SSA и ещё больше с приложениями SSA. Введение в SSA может быть найдено в Elsner and Tsonis (1996). Более углубленными работами являются монография Golyandina et al. (2001) (её содержание кратко и частично изложено на русском языке в учебных пособиях (Голяндина, 2004)), обзор Ghil et al. (2002), специальный выпуск журнала ‘Statistics and Its Interface’ (Zhigljavsky, 2010, Ed.) и книга Golyandina and Zhigljavsky (2013). В книге Golyandina et al. (2018) описаны и структурированы различные обобщения и модификации метода SSA для анализа временных рядов и изображений, а также их алгоритмы и их реализация в R-пакете Rssa. Сайт-компаньон к книге с примерами на R — http://ssa-with-r-book.github.io

Литература

• Голяндина Н. Э. (2004) Метод «Гусеница»-SSA: анализ временных рядов: Учеб. пособие. СПб: Изд-во СПбГУ.
• Голяндина Н. Э. (2004) Метод «Гусеница»-SSA: прогноз временных рядов: Учеб. пособие. СПб: Изд-во СПбГУ.
• Данилов Д. и Жиглявский А. (Ред.) (1997): Главные компоненты временных рядов: метод «Гусеница», Изд. СПбГУ.
• http://www.gistatgroup.com/gus/papers.html
• Ефимов В. М., Галактионов Ю. К. О возможности прогнозирования циклических изменений численности млекопитающих // Ж. общ. биол., 1983. № 3, с.343-352.
• Ефимов В. М., Галактионов Ю. К., Шушпанова Н. Ф. Анализ и прогноз временных рядов методом главных компонент. -Новосибирск: Наука, 1988. −70с.

Англоязычная литература

• Broomhead, D.S., and G.P. King (1986a): "Extracting qualitative dynamics from experimental data", Physica D, 20, 217—236.
• Broomhead, D.S., and G. P. King (1986b): "On the qualitative analysis of experimental dynamical systems". Nonlinear Phenomena and Chaos, Sarkar S (Ed.), Adam Hilger, Bristol, 113—144.
• Elsner, J.B. and Tsonis, A.A. (1996): Singular Spectrum Analysis. A New Tool in Time Series Analysis, Plenum Press.
• Fraedrich, K. (1986) "Estimating dimensions of weather and climate attractors". J. Atmos. Sci. 43, 419—432.
• Ghil, M., and R. Vautard (1991): "Interdecadal oscillations and the warming trend in global temperature time series", Nature, 350, 324—327.
• Ghil, M., R. M. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, et al. (2002) "Advanced spectral methods for climatic time series", Rev. Geophys. 40(1), 3.1-3.41.
• Golyandina, N., A. Korobeynikov and A. Zhigljavsky (2018): Singular Spectrum Analysis with R. Series Use R! Springer Verlag. ISBN 3662573784.
• Golyandina, N., V. Nekrutkin and A. Zhigljavsky (2001): Analysis of Time Series Structure: SSA and related techniques. Chapman and Hall/CRC. ISBN 1-58488-194-1.
• Golyandina, N., and E. Osipov (2007) "The ‘Caterpillar’-SSA method for analysis of time series with missing values", J. Stat. Plan. Inference 137(8), 2642—2653.
• Golyandina, N. and D. Stepanov (2005): "SSA-based approaches to analysis and forecast of multidimensional time series". In: Proceedings of the 5th St.Petersburg Workshop on Simulation, June 26-July 2, 2005, St. Petersburg State University, St. Petersburg, pp. 293–298.
• Golyandina, N. and K. Usevich (2010): "2D-extension of Singular Spectrum Analysis: algorithm and elements of theory". In: Matrix Methods: Theory, Algorithms and Applications (Eds. V.Olshevsky and E.Tyrtyshnikov). World Scientific Publishing, 449—473.
• Golyandina, N., and A. Zhigljavsky (2013) Singular Spectrum Analysis for time series. Springer Briefs in Statistics, Springer, ISBN 978-3-642-34912-6.
• Kondrashov, D., and M. Ghil (2006): "Spatio-temporal filling of missing points in geophysical data sets", Nonlin. Processes Geophys., 13, 151—159.
• Nekrutkin, V. (2010) "Perturbation expansions of signal subspaces for long signals". J. Stat. Interface 3, 297—319.
• de Prony, G. (1795) "Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques et sur celles de la force expansive de la vapeur de l’eau et la vapeur de l’alkool à différentes températures". J. de l’Ecole Polytechnique, 1(2), 24-76.
• Vautard, R., and M. Ghil (1989): "Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series", Physica D, 35, 395—424.
• Zhigljavsky, A. (ed.) (2010) Statistics and Its Interface (Special issue on the singular spectrum analysis in time series), vol 3. Guest Editor.

См. также

• Метод главных компонент

Ссылки

• «Гусеница»-SSA Статьи и программы от Gistat Group.
• Rssa Эффективная реализация SSA в R. Сайт с примерами.
• Бесплатная программа SSA-MTM Toolkit для Unix от UCLA.
• kSpectra Toolkit для Mac OS X от SpectraWorks.
• Ещё одна SSAwiki страница.