close

Вход

Забыли?

вход по аккаунту

?

Dik

код для вставкиСкачать
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное автономное
образовательное учреждение высшего образования
САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
АЭРОКОСМИЧЕСКОГО ПРИБОРОСТРОЕНИЯ
О. Е. Дик
НЕЛИНЕЙНЫЙ АНАЛИЗ
ВРЕМЕННЫХ РЯДОВ
Учебное пособие
УДК 51-77
ББК 22.17
Д45
Рецензенты:
доктор физико-математических наук, профессор В. Г. Фарафонов;
доктор физико-математических наук, профессор Ю. А. Пичугин
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
Дик, О. Е.
Д45 Нелинейный анализ временных рядов: учеб. пособие /
О. Е. Дик. – СПб.: ГУАП, 2018. – 51 с.
ISBN 978-5-8088-1254-3
Включает описание необходимых теоретических сведений и современных методов анализа временных рядов (корреляционного,
вейвлетного, мультифрактального и рекуррентного), позволяющих
оценивать динамические изменения в паттернах сложных нестационарных сигналах, порождаемых физическими, химическими, биологическими или экономическими системами.
Составлено в соответствии с программой спецкурса «Нелинейный
анализ временных рядов» для студентов-магистров очного отделения
факультета инноватики и базовой магистерской подготовки СанктПетербургского государственного университета аэрокосмического
приборостроения.
УДК 51-77
ББК 22.17
ISBN 978-5-8088-1254-3
©
©
Дик О. Е., 2018
Санкт-Петербургский государственный
университет аэрокосмического
приборостроения, 2018
1. ПОНЯТИЕ ВРЕМЕННОГО РЯДА,
ДИНАМИЧЕСКОЙ СИСТЕМЫ
И АТТРАКТОРА ДИНАМИЧЕСКОЙ СИСТЕМЫ
Изучение сложных физических, химических, биологических или
экономических систем часто реализуется посредством изучения сигналов, порождаемых этими системами. Так, для анализа состояния
финансовых рынков используются данные о колебаниях цен и объемов продаж ценных бумаг, для анализа функционального состояния человека – данные его электроэнцефалограммы и электрокардиограммы, анализа состояния земной коры – данные о ее колебаниях, полученные в различные промежутки времени. Все эти сигналы
представляют собой скалярные временные ряды, которые, в свою
очередь, являются последовательностями значений некоторых переменных, регистрируемых через определенные интервалы времени:
{x(ti)}, ti = t0 + (i – 1)dt, i = 1, …, n.
При анализе временного ряда решают две основные задачи: задачу идентификации и задачу прогноза [Лоскутов, Михайлов, 2007].
Решение задачи идентификации означает получение ответа на вопрос: каковы параметры системы, породившей изучаемый временной
ряд (то есть какие величины размерности вложения, энтропии, корреляционной и фрактальной размерностей характеризуют данный ряд,
каковы его рекуррентные свойства) [Лоскутов, Михайлов, 2007].
Решение задачи прогноза означает получение возможности предсказания по данному временному ряду его значений в будущие промежутки времени (это может быть, например, реализовано на основании выявления долговременных корреляций в исследуемом временном ряду) [Лоскутов, Михайлов, 2007].
1.1. Свойства временных рядов: стационарность,
нестационарность, стохастичность, детерминизм
Большинство изучаемых биологических, химических и экономических систем проявляют сложную нелинейную динамику. Это
означает, что их колебания сильно вариабельны, а связи между
причинами этой вариабельности и вызванными ими последствиями существенно нелинейны и, соответственно, малые нерегулярные
флуктуации могут вызвать значительные изменения в паттернах
колебаний [Ivanov, et al., 1999].
Стационарность – свойство системы не менять свои характеристики со временем. Соответственно, нестационарные временные
3
ряды изменяют со временем свои характеристики (корреляционную и фрактальную размерности, рекуррентные свойства и т. д.).
Стохастичность (случайность) – свойство системы не быть детерминированной. В случае стохастического временного ряда последующее состояние системы описывается как величинами, которые
могут быть предсказаны, так и случайными.
1.2. Понятие аттрактора динамической системы
Аттрактор – это множество состояний (точек фазового пространства) динамической системы, к которому она стремится с течением времени. Другими словами, динамическая система, эволюционируя, приходит в определенную область фазового пространства, а все траектории
системы стягиваются к ее аттрактору [Безручко, Смирнов, 2005].
Наиболее простыми вариантами аттрактора являются притягивающая неподвижная точка и периодическая траектория (предельный цикл). Большинство типов движения описывается простыми
аттракторами, являющимися ограниченными циклами. Например, фазовый портрет периодического процесса – предельный цикл,
и размерность этого аттрактора целочисленна.
Хаотическое движение описывается странными аттракторами, которые очень сложны и имеют много параметров. Например, простая
трехмерная система погоды описывается известным аттрактором Лоренца. Размерность странного аттрактора не- целочисленна (меньше
размерности его фазового пространства), то есть фрактальна. Таким
образом, среди аттракторов можно выделить неподвижные точки,
предельные циклы и хаотические аттракторы [Kantz, Schreiber, 1997].
Если поведение системы описывается системой дифференциальных уравнений, то размерность фазового пространства равна числу
переменных системы, а размерность аттрактора меньше или равна
числу этих переменных.
Чем больше переменных описывает поведение динамической системы, тем более сложным является это поведение, и, соответственно, менее предсказуемо. Размерность аттрактора в этом случае также больше.
1.3. Реконструкция аттрактора из исходного временного ряда
Одним из методов визуализации динамического поведения
сложных нелинейных систем является представление их поведения в фазовом пространстве. Это позволяет обнаружить изменения
во времени значений, которые принимают переменные системы.
4
Визуализация и связанная с ней реконструкция аттрактора
динамической системы, основана на теореме вложения Такенса
[Takens, 1981], согласно которой фазовая траектория n – мерной динамической системы лежит на аттракторе, принадлежащем гладкому k – мерному многообразию. Тогда по одномерной проекции x(t)
этой траектории методом задержки можно получить m – мерную
реконструкцию исходного аттрактора как множество векторов
y(t) = (x(t), x(t+d), …, x(t + (m – 1)d), m ≥ 2k + 1,
где d – временная задержка, m – размерность вложения (минимальная размерность пространства, в котором восстановленный аттрактор воспроизводит свойства исходного аттрактора, то есть наименьшее число независимых переменных, однозначно определяющих
его свойства).
Таким образом, каждому значению x(ti) временного ряда {x(t)} в заданный момент времени соответствует точка (x(ti), x(ti + d), x(ti + 2d),
x(ti + (m – 1)d)), которая определяет значения той же переменной после некоторого фиксированного времени задержки. Множество таких
точек для последовательных моментов времени образует траекторию,
описывающую эволюцию исследуемой системы (временного ряда)
[Kantz, Schreiber, 1997].
Следовательно, фазовый портрет динамичекой системы может
быть построен по временному ряду, представляющему собой зависимость от времени только одной из наблюдаемых переменных
системы (а не всех ее переменных), если в качестве недостающих
координат вектора состояний использовать исходный ряд, взятый
с некоторым запаздыванием.
В соответствии с теоремой Такенса, можно подобрать такие значения параметров m и d, что в результате преобразования множество полученных точек будет воспроизводить аттрактор исследуемой динамической системы [Kantz, Schreiber, 1997].
Пространство, задаваемое для восстановления аттрактора, называется пространством вложения или лаговым пространством. Размерность восстановленного аттрактора является мерой стохастичности исследуемого процесса, так как чем меньше эта размерность,
тем сильнее детерминирован процесс, а чем больше размерность,
тем более стохастичен процесс.
Для корректной реконструкции фазовой траектории состояний
y(t) из исходного временного ряда необходимо выбрать такую задержку d, чтобы координаты точки (x(ti), x(ti + d)) несли как можно
меньше повторяющейся информации.
5
Подбор оптимального интервала задержки производится на основании поиска первого минимума функции взаимной информации [Fraser, Swinney, 1986]. Далее методом поиска ближайших
ложных соседей определяется оптимальная размерность вложения m [Kennel et al.,1992].
Функция взаимной информации вычисляется следующим образом:
S(d) = −∑ pij (d)ln
pij (d)
i,j
pi pj
,
где pi – вероятность, с которой элемент временного ряда может оказаться в i-м интервале, pj- вероятность, с которой элемент временного ряда может оказаться в j-м интервале при разбиении промежутка, равного размаху амплитуды сигнала, на несколько интервалов
[Fraser, Swinney, 1986]. В этом случае pij(d) – совместная вероятность того, что один элемент временного ряда окажется в i-м интервале, а другой, взятый с задержкой d в j-м интервале.
Алгоритм поиска ложных ближайших соседей состоит в следующем:
1) восстанавливается аттрактор в лаговом пространстве и для
каждой его точки осуществляется поиск ее ближайших соседей;
2) далее аттрактор восстанавливается в лаговом пространстве,
размерность которого на единицу больше и определяется число
ложных ближайших соседей. Ложными ближайшими соседями будут пары точек, для которых отношение между расстояниями
R=
xi +1 − xj +1
xi − xj
до и после увеличения размерности лагового пространства (размерности вложения m) превышает некоторый порог.
Отметим, что при увеличении значения m число ложных соседей
уменьшается, а при правильном выборе размерности вложения число ложных соседей будет минимальным (равным нулю).
Величина размера окрестности ε точки x выбирается равной 1%
от величины стандартного отклонения анализируемого временного
ряда.
На рис. 1 показан сигнал, несущий две частоты. Автокорреляционная функция для этого сигнала имеет первый минимум при
значении временной задержки d = 25, а функция взаимной информации имеет первый минимум при значении временной задержки
d = 15. В связи с этим для двух значений параметра d далее опре6
d = 15
0
–1
10
8
6
4
2
t (с)
Автокорреляционная функция
1
0,5
0
–0,5
0
10
20
d
0,5
0
1
1,5
2
2,5
3
d = 25
0,5
0
40
30
Фракция ложных соседей
Двухчастотный сигнал
1
1
3
1,5
2
2,5
Размерность вложения m
Функция взаимной информации
2,5
2
1,5
1
0,5
0
5
10
d
15
20
Рис. 1. Определение размерности вложения m и временной задержки d
для двухчастотного сигнала
делим оптимальную размерность вложения m методом поиска ближайших ложных соседей. Как следует из рис. 1, величина m = 2 для
обоих значений параметра d.
Далее построим портреты фазовых траекторий в реконструированном фазовом пространстве для анализируемого двухчастотного
сигнала для размерности вложения m = 2 и различных значениях
параметра d (Рис. 2).
Как видно из рис. 2, оптимальным является величина задержки
d = 15.
Портреты фазовых траекторий в реконструированном фазовом
пространстве для зашумленного периодического сигнала с периодом 1 и стохастического сигнала, описывающего движение броуновской частицы, а также для колебаний зашумленной системы
Лоренца, проявляющей зашумленное хаотическое поведение, представлены на рис. 3. Очевидно, что фазовые траектории даже зашумленных периодического и хаотического сигналов проявляют хорошо выраженную структуру, в отличие от траектории, полученной
для стохастического сигнала.
7
x(t + d)
m = 2, d = 25
m = 2, d = 15
1
1
0,5
0,5
0
0
–0,5
–0,5
–1
–1
–1 –0,5
0
0,5
–1 –0,5
1
m = 2, d = 10
0 0,5
x(t)
1
x(t + d)
1
0,5
0
–0,5
–1
–1 –0,5
0 0,5
x(t)
1
Рис. 2. Портреты фазовых траекторий в реконструированном фазовом
пространстве для двухчастотного сигнала, построенные
при различных значениях параметров d и m
x(t + 2d)
x(t + 2d)
2
0
–2
2
x(t + 2d)
Зашумленный хаотический
20
0
–20
20
2
0
–2
2
10
0
–10
–20 –20
–10
0
10
20
0
1
2
0
1
2
Зашумленный периодический
1
0
–1
–2 –2
–1
Стохастический
1
0
x(t + d)
–1
–2 –2
–1
x(t)
Рис. 3. Портреты фазовых траекторий в реконструированном фазовом
пространстве для зашумленных хаотического, периодического
и стохастического сигналов. Величины параметров d = 5 и m = 3
для хаотического, d = 15 и m = 3 для периодического , d = 40 и m = 9
для стохастического сигналов
8
Таким образом, структура реконструированного аттрактора хорошо видна только при правильном выборе параметров реконструкции аттрактора m и d. При этом, если d слишком мало, то i и i + 1 координаты точки в фазовом пространстве практически неотличимы.
Если d слишком велико, то координаты оказываются некореллированными. И в том, и в другом случае реконструированный аттрактор не будет отражать истинной динамики процесса.
Вопросы:
1. Что такое временной ряд и каковы его свойства?
2. Что такое аттрактор динамической системы и какие типы аттракторов существуют?
3. В чем состоит реконструкция аттрактора динамической системы?
4. Как подбираются параметры реконструкции аттрактора для
заданного временного ряда?
5. Как найти оптимальную величину лага (интервала задержки)?
6. Зачем нужна автокорреляционная функция и функция взаимной информации при построении аттрактора?
7. Как найти оптимальную размерность вложения?
Лабораторная работа 1. Задать систему дифференциальных
уравнений (например, систему уравнений Лоренца
dx
=
σ(y − x)
dt
dy
= rx − y − xz
dt
dz
=
−bz + xy
dt
с параметрами σ = 10, r = 28, b = 8/3), найти решение этой системы
и построить график. Изменить значения параметра и построить зависимость решения от значения параметра.
Лабораторная работа 2. Написать программу, моделирующую
хаотический аттрактор, которым обладает отображение Хенона
xn +1 =
1 − axn2 + byn
yn +1 = xn
с параметрами a = 1,4, b = 0,3 и построить ее решение.
Лабораторная работа 3. Написать программу реконструкции
аттрактора по заданному временному ряду, найти оптимальную ве9
личину лага и размерности вложения. Визуализировать реконструированный аттрактор.
Лабораторная работа 4. Написать программу, моделирующую
хаотический аттрактор, которым обладает система дифференциальных уравнений Рёсслера
dx
=−y − z
dt
dy
= x + ay
dt
dz
=b − cz + xy
dt
с параметрами a = 0,2, b = 0,2, c = 5 и построить ее решение.
10
2. ВЫЧИСЛЕНИЕ КОРРЕЛЯЦИОННОГО ИНТЕГРАЛА
И КОРРЕЛЯЦИОННОЙ РАЗМЕРНОСТИ
ВОССТАНОВЛЕННОГО АТТРАКТОРА
Корреляционным интегралом называется усредненная вероятность того, что два состояния системы x(i) и x(j) в два различных
момента времени окажутся близкими. Оценка корреляционного интеграла дается корреляционной суммой:
=
C(ε) lim
1
n
∑ Θ(ε −
2
n →∞ n i,j =1
)
x=
i − xj , i, j 1, ..., n,
где n – количество рассматриваемых состояний; ε – размер окрестности точки x, ||...|| – норма, q (·) – функция Хевисайда (принимающая значения 0 или 1) [Kantz, Schreiber, 1997].
Фактически, корреляционная сумма определяет число точек, отстоящих от данной точки на расстояние, не превышающее ε.
Корреляционный интеграл используется для вычисления корреляционной размерности, которая является количественной характеристикой аттрактора, несущей информацию о степени сложности
поведения динамической системы,
2.1. Алгоритм вычисления
корреляционной размерности аттрактора
Если зависимость C(e) имеет степенной вид
C(ε) ≈ ε D ,
тогда исследуемый временной ряд обладает свойством фрактальности, а величина D
называется его корреляционной размерностью.
Из приведенной выше формулы следует, что корреляционная
размерность восстановленного аттрактора мoжет быть вычислена
следующим образом
ln C(ε)
D = lim
.
ε→0 ln ε
Алгоритм вычисления корреляционной размерности аттрактора
состоит из следующей последовательности процедур:
1) вычисляется корреляционная сумма, 2) выделяется на графи)) D ln(ε) область линейной зависимости и находится знаке ln(C(ε=
чение тангенса угла наклона прямой.
11
Отметим, что поиск линейного участка зависимости корреляционного интеграла от величины измельчения разбиения (ε) в двойном логарифмическом масштабе может быть затруднен при анализе
сложного или сильно зашумленного сигнала. В этом случае используют оценку Такенса – Тейлера [Kantz, Schreiber, 1997]:
DT (ε) =ε
C(ε)
.
∫ (C(x)) / x)dx
0
Поскольку восстановление аттрактора можно рассматривать как
проецирование во вложенное (лаговое) пространство, то величина D
должна быть меньше m. Если при этом продолжать увеличивать m,
то начиная с некоторого значения m, должно наступить насыщение
(величина D должна перестать увеличиваться). Это является признаком детерминистского процесса, если же исследуемый процесс
является стохастическим, то насыщение на графике D(e) при увеличении m не происходит [Kantz, Schreiber, 1997].
Кроме того, необходимо отметить следующее. Во-первых, для правильного вычисления корреляционной размерности длина (n) временного ряда не может быть меньше чем n > 10D/2. Во-вторых, точки, расположенные в ряду в непосредственной близости друг от друга, являются сильно коррелированными. Для исключения таких корреляций и
уменьшения ошибок в вычислении корреляционного интеграла не следует принимать в расчет точки, расположенные в последовательности
на расстоянии, меньшем, чем величина w, называемая окном Тейлера.
Размер окна Тейлера определяется из графика пространственновременного разделения, то есть кривой, описывающей плотность
вероятности того, что две точки временного ряда, находящиеся на
расстоянии dt, окажутся в восстановленном аттракторе на расстоянии, меньшем величины e [Kantz, Schreiber, 1997].
При этом, как правило, строится семейство кривых пространственновременного разделения, соответствующих разным значениям плотности
вероятности, а размер окна Тейлера определяется как величина, соответствующая первому локальному минимуму, общему для всех кривых.
На рис. 4. показан хаотический сигнал, соответствующий решению системы дифференциальных уравнений Лоренца, а также
графики зависимостей D(e) при различных значениях параметра m.
Параметры аттрактора этой системы получены на основе вычислений минимумов автокорреляционной функции, функции взаимной
информации и фракции (относительного числа) ложных соседей,
12
Cигнал
Фракция ложных соседей
20
0,5
0
–20
0
100
tau
0
300
200
Автокорреляционная функция
0
–1
0
100
50
d
5
2
3
4
Hазмерность вложения m
Оцениватель Такенса-Тейлера
10
D
1
1
5
0
10–2
100
Расстояние e
102
Функция взаимной информации
2,5
2
1,5
1
0,5
5
0
10
d
15
20
Рис. 4. Вычисление корреляционной размерности ряда,
соответствующего решению системы уравнений Лоренца
Видно, что графики имеют четко выраженное плато в широком диапазоне значений e для всех значений размерности вложения m. Это
позволяет определить корреляционную размерность анализируемого временного ряда D = 2,1 (для m > 1).
Вопросы:
1. Как вычислить корреляционную размерность аттрактора?
2. В чем состоят ограничения на величину e при вычисления корреляционного интеграла?
3.В чем состоит сложность вычисления корреляционной размерности аттрактора?
Лабораторная работа 1. Вычислить корреляционный интеграл
и найти величину корреляционной размерности аттрактора для
отображения Хенона
xn +1 =
1 − axn2 + byn
yn +1 = xn
с параметрами a = 1,4, b = 0,3.
13
Лабораторная работа 2. Вычислить корреляционную размерность аттрактора для системы уравнений Рёсслера
dx
=−y − z
dt
dy
= x + ay
dt
dz
=b − cz + xy
dt
с параметрами a = 0,2, b = 0,2, c = 5 и выявить влияние параметров
на точность вычисления для заданного временного ряда.
14
3. ВЕЙВЛЕТ-АНАЛИЗ ВРЕМЕННЫХ РЯДОВ
Вейвлет-преобразование широко используется для решения
широкого класса задач, связанных с анализом временных рядов,
например, анализом биологических сигналов (электрокардиограмм, электроэнцефаллограмм и томограмм мозга) [Астафьева,
1996, Дремин и др., 2001], а также анализом финансовых временных рядов (изменением цен и объемов продаж ценных бумаг) [Ширяев, 2007].
Метод вейвлет-преобразования [Percival, 2000] позволяет разложить одномерный сигнал, представленный в виде ряда числовых
N
значений {x(ti )}i =1 , по набору копий одной исходной функции – прототипа (материнского вейвлета ψ(t)). Эти копии (вейвлетные функции) масштабируются в некоторое число раз (то есть растягиваются
во времени) и смещаются во времени на некоторое расстояние:
ψ a,t0 =
1  t − t0 
ψ
,
a  a 
где a и t0 – параметры масштаба и сдвига, определяющие ширину
и смещение вейвлета. При этом временной масштаб (a) обратно пропорционален частоте преобразования Фурье (f).
Смещение вейвлета вдоль изучаемого сигнала дает возможность
обнаружить изменение масштаба, а значит и частоты сигнала во
времени, что важно для существенно нестационарных рядов, какими являются большинство анализируемых физических или экономических рядов [Малла, 2005].
Подобно Фурье – преобразованию, вейвлет – преобразование непрерывного сигнала представляет собой свертку сигнала с вейвлетной функцией:
+∞
W (=
a,t0 )
*
∫ x(t)ψa,t0 dt,
−∞
ψ*a,t0
где символ
означает комплексно-сопряженную функцию. Величина коэффициента W(a, t0) вейвлет – преобразования характеризует наличие и интенсивность соответствующего временного масштаба a в момент времени t0.
3.1. Дискретное вейвлет преобразование временного ряда
Для дискретного сигнала вейвлет преобразование вычисляется для дискретных значений параметров масштаба a = 2m и сдви15
га t0 = k2m, где k, m – целые числа. В этом случае семейство вейвлетных функций имеет следующий вид [Добеши, 2001]:
ψm,k (t) =
 t − k2m 
ψ
.
m 

2m  2

1
При этом дискретное преобразование сигнала на m – уровне разложения представляет собой суперпозицию вейвлетных функций.
После обратного вейвлет-преобразования анализируемый сигнал x(t) равен сумме сглаженной (низкочастотной) компоненты последнего уровня (Am) и детализирующих (высокочастотных) компонент всех уровней (Dm, …, D1) разложения (рис. 5):
x(=
t) Am (t) + Dm (t) + ... + D1 (t).
Компонента A3(t) описывает грубое приближение к исходному сигналу на третьем уровне разложения, а D1(t), … D3(t) – мелкомасштабные детали, полученные на трех уровнях разложения; причем компонента D1(t)
характеризует наиболее высокочастотные детали анализируемого ряда.
Таким образом, дискретное вейвлет преобразование фактически
является прохождением сигнала через определенные частотные
фильтры, а кратномасштабный алгоритм многоуровнего разложения и последующего восстановления сигнала позволяет выявить
особенности этого сигнала на различных частотах [Малла, 2005].
Пример сглаженной и детализирующих компонент вейвлет-разложения сигнала представлен на рис. 6.
В отличие от оконного преобразования Фурье, имеющего одно и
то же разрешение по времени и частоте для всех точек сигнала, вейвлетные функции имеют уменьшающееся (при уменьшении масштаба) разрешение по времени и увеличивающееся (при увеличении масштаба) разрешение по частоте. Это свойство вейвлет – преобразования дает ему большое преимущество при анализе сигнала, так как
x(t)
A1(t)
A2(t)
A3(t)
D1(t)
D2(t)
D3(t)
Рис. 5. Трехуровневое вейвлет-разложения сигнала:
x(t) = A3 (t) + D3 (t) + D2 (t) + D1 (t)
16
Cигнал
0,5
0
–0,5
200
400
600
800
1000
0
200
400
600
800
1000
0
200
400
600
800
1000
0
200
400
600
800
1000
0
200
400
600
800
1000
A3
0,5
0
–0,5
0
D3
0,5
0
–0,5
D1
D2
0,5
0
–0,5
0,5
0
–0,5
nL
Рис. 6. Трехкомпонентное вейвлет-разложение сигнала. Сверху вниз
показаны сигнал и сглаженная A3, а также детализирующие
компоненты D1, D2, D3
«адаптивное» частотно – временное окно позволяет с хорошей точностью извлекать высокочастотную информацию из относительно малых временных интервалов, а низкочастотную информацию определять из относительно широких временных интервалов сигнала.
Для анализа полученных компонент сигнала обычно используется
быстрый алгоритм дискретного преобразования Фурье с последующим оцениванием спектральной плотности энергии этих компонент,
что позволяет определять распределение энергии компонент сигнала
или их мощности по частотам [Смоленцев, 2005; Блаттер, 2006].
3.2. Интегральная оценка
спектральной плотности энергии временного ряда
Спектральная плотность S(f) энергии компоненты D(t) временного ряда равна квадрату Фурье преобразования этой компоненты и
описывает распределение энергии по частоте:
2
S(f ) = ∫ D(t)e−2πiftdt .
17
Спектральная плотность энергии
0,16
Интегральная энергия
2
1,8
0,14
fmax
1,6
0,12
1,4
0,1
1,2
0,08
1
0,06
0,8
0,6
0,04
f1
0,02
0
Emax
0,4
f1
0,2
0
5
частота (Гц)
0
10
0
5
10
частота (Гц)
15
Рис. 7. Определение параметров спектральной плотности
энергии компоненты сигнала
Интеграл от cпектральной плотности энергии компоненты D(t)
f2
E = ∫ S(f )df
f1
определяет накопление спектральной плотности энергии данной
компоненты в пределах заданной полосы частот [f1, f2].
В качестве параметра оценки спектральной плотности в работе [Мусалимов и др., 2009] предложена следующая величина:
Emax
k=
,
(f2 − f1 )fmax
где fmax – значение частоты, соответствующее максимуму спектральной плотности, Emax – максимальный уровень накопления спектральной плотности энергии, [f1, f2] – интервал значений частот, оставшийся после 5% фильтрации шумовых компонент, то есть f1 и f2 соответствуют значениям 0,05*Smax и 0,95*Smax (рис. 7).
В приведенной формуле параметр k определяет связь между
уровнем максимального накопления энергии компоненты сигнала, частотой, соответствующей максимальному накоплению спектральной плотности энергии, и интервалом частот, на котором накоплена эта энергия.
18
3.3. Непрерывное вейвлет преобразование временного ряда
и оценивание энергии сигнала
Непрерывное вейвлет преобразование временного ряда, описывающего исследуемый сигнал x(t), обычно записывают в виде:
+∞
=
W (a,t0 )
1
 t − t0 
x(t) ψ* 
 dt,
∫
a −∞
 a 
где a – параметр масштаба, t0 – параметр временного сдвига,
ψ((t – t0)/a) – вейвлет функция, полученная из базисного вейвлета
ψ(t) путем масштабирования и сдвига по времени, символ * означает
комплексное сопряжение [Добеши, 2001].
3.4. Базисный вейвлет Морле
В качестве базисного вейвлета при исследовании сигналов, как правило, применяют комплексный вейвлет Морле [Grossmann et al., 1984]:
(
)
ψ0 (t) =
π−1/4 exp(−0,5t2 ) exp(iw0t) − exp(−0,5w20 ) ,
где второй компонентой в скобках при w0 = 2p > 0можно пренебречь.
Значение ω0 = 2π обеспечивает простое соотношение между масштабом a и частотой f спектра Фурье: f = 1/a.
Тогда вейвлет преобразование сигнала x(t) можно записать в виде:
W (f,t0 ) =π−1/4 f
+∞
2 2
∫ x(t)exp(−0.5(t − t0 )
−∞
f )exp(−i2π(t − t0 )f )dt.
Величина модуля вейвлетного спектра |W(f, t0)| характеризует
наличие и интенсивность частоты f в момент времени t0 в анализируемом сигнале, а величина квадрата модуля |W(f, t0)|2 описывает
мгновенное распределение энергии фрагмента сигнала по частотам,
то есть локальный вейвлетный спектр энергии в момент времени t0.
Этот спектр характеризует локализацию определенной частоты по
времени ее возникновения в исследуемом ряду.
3.5. Глобальный вейвлетный спектр.
После интегрирования локальных спектров по времени можно
получить глобальный вейвлетный спектр, то есть
E(f ) =
t2
∫ W (f,t0 ) dt0
2
t1
19
Cигнал
x
50
0
–50
0
1
2
3
4
5
t (c)
6
7
8
9
7
8
9
Локальный вейвлетный спектр
f(Гц)
15
10
5
E(f)
10
0
x106
1
2
3
4
5
t0 (с)
6
Глобальный вейвлетный спектр
5
0
8
f(Гц)
16
Рис. 8. Примеры локального и глобального спектров сигнала
интегральное распределение энергии вейвлетного спектра сигнала
по частотам в определенном интервале времени [t1, t2].
Примеры локального и глобального спектров сигнала приведены
на рис. 8.
На представленном рисунке глобальный вейвлетный спектр показывает широкий интервал частот в исследуемом сигнале с максимумом на частоте 11 Гц. Однако локальный вейвлетный спектр демонстрирует, что в данном интервале записи максимальная частота
сигнала присутствовала не постоянно, а в некоторые определенные
и выделяемые в этом локальном спектре промежутки времени. Отметим, что спектр мощности Фурье информацию о временной локализации максимальной частоты в сигнале не дает, как не дает ее и
его аналог (глобальный вейвлетный спектр).
На рис. 9 изображен периодический сигнал x(t), частота которого
первые 2 сек составляла 3,84 Гц, а последующие – увеличилась до
62,8 Гц, то есть заданный следующей функцией:
 sin(2πw
=
1t), w1 3, åñëè 0 < t < 2,
x(t) = 
w2t), w2 10, åñëè 2 ≤ t < 4,
 sin(2π=
20
Сигнал
x
2
0
–2
0
0,5
1,5
1
2
t (c)
2,5
3
3,5
3
3,5
Локальный вейвлетный спектр
f (Гц)
20
10
0
0
0,5
1,5
1
2
t0 (с)
2,5
E (f)
Глобальный вейвлетный спектр
150
100
50
0
0
10
f(Гц)
5
20
15
Рис. 9. Вейвлетные спектры сигнала
Сигнал
x
2
0
–2
0
0.5
1
1.5
2
t (c)
2.5
3
3.5
3
3,5
Локальный вейвлетный спектр
f (Гц)
20
10
0
0
0,5
1
1,5
2
t0 (с)
2,5
Глобальный вейвлетный спектр
E (f)
50
0
0
2
4
6
f (Гц)
8
10
12
Рис. 10. Вейвлетные спектры сигнала
21
В отличие от этого, на рис. 10 показан двухпериодический сигнал y(t), содержащий обе частоты 3,84 Гц и 62,8 Гц одновременно,
то есть заданный функцией:
y(t) = sin(2πw1t) + 0.5 sin(2πw2t), w1 = 3, w2 = 10.
Из приведенных рисунков видно, что локальный вейвлетный
спектр демонстрирует различия в частотно-временных особенностях сигналов x(t) и y(t).
3.6. Удаление шумовой составляющей сигнала
Метод вейвлет преобразования применяется также для удаления шумовой составляющей из анализируемого сигнала (так называемый wavelet denoising method) [Малла, 2005].
В случае использования вейвлет-преобразования пороговая обработка вейвлет коэффициентов осуществляется по принципу
Штейна несмещенной оценки риска. При этом используется схема
Регистрируемый сигнал и тренд (denoising method)
1
0
–1
0
5
10
15
20
25
30
20
25
30
(detrended method)
1
0
–1
0
5
10
15
t (c)
Рис. 11. Сравнение эффективности удаления тренда
из оригинальной записи данных с помощью адаптивного метода
удаления тренда и с использованием вейвлетов. Регистрируемый сигнал
показан тонкой линией, тренд – жирной линией
22
мягкой фильтрации, при которой коэффициенты дискретного вейвлет-преобразования dj,k изменяются по формуле:
0

dj,k = 
sign (d j,k ) ⋅

(
d j,k ≤ p
d
j,k
−p
)
dj,k > p
где p – выбранное значение порога на масштабе j связано с числом N
точек оцифровки исходной функции и дисперсией вейвлет-коэффициентов σ на первом масштабе следующим образом
σ = 2∑ d12,k / N [Малла, 2005].
k
Альтернативным подходом является адаптивный метод удаления тренда (the adaptive detrending method) [Hu, et al, 2009]. Адаптивный метод удаления тренда основан на разделении временного
ряда на сегменты длины w с последующей аппроксимацией каждого сегмента полиномом порядка m.
Рис. 11 иллюстрирует колебательные тренды, выделенные из
исследуемого временного ряда с использованием двух описанных
выше алгоритмов.
На основании сравнения полученных трендов можно заключить,
что адаптивный метод удаления тренда более аккуратно удаляет
медленный колебательный тренд.
3.7. Сравнение вейвлетного преобразования сигнала
с быстрым преобразованием Фурье
и оконным преобразованием Фурье
Преимущества вейвлетного анализа (CWT) сигналов, по сравнению с классическим спектральным анализом, основанном на дискретном преобразовании Фурье, широко известны.
В случае, когда необходимо только определить усредненные по
времени спектральные компоненты, представленные в анализируемом сигнале, классический подход может быть успешно применен.
Однако если нас интересует распределение во времени различных
частотных компонент, использование вейвлетов имеет значительные
преимущества. В отличие от быстрого преобразования Фурье (FFT),
непрерывное вейвлет преобразование позволяет получить информацию об изменении частотных характеристик сигнала во времени,
то есть узнать, в какие именно моменты времени возникают те или
иные частоты [Малла, 2005].
23
Другим подходом к исследованию сигналов может быть использование оконного преобразования Фурье (STFT) [Allen, 1977].
Алгоритм этого метода состоит из следующей последовательности процедур:
1) регистирируемый сигнал разделяется на короткие перекрывающиеся фрагменты,
2) для каждого фрагмента применяется быстрое преобразование
Фурье с оконной функцией, например, окном Хэмминга,
3) далее окно сдвигается на число значений, равных разности w –
n, где w – длина окна, n – число значений, которые перекрываются
в каждом сегменте.
В силу такого алгоритма спектр, полученный при помощи оконного преобразования Фурье, является произведением сигнала и
оконной функции, а результат оконного преобразования Фурье зависит от обоих параметров w и n. Поэтому при использовании оконного преобразования Фурье оказывается невозможным одновременно обеспечить хорошее разрешение как по времени, так и по частоте. Чем уже окно, тем выше разрешение по времени, но ниже по
частоте, то есть узкое окно не обеспечивает хорошего определения
частот в структуре исследуемого сигнала. Таким образом, фиксированный размер окна в методе STFT не позволяет описать локальные свойства паттернов сигнала, так как увеличение разрешения по
времени приводит к ухудшению разрешения по частоте.
В отличие от оконного преобразования Фурье, в вейвлетном анализе используются окна различных размеров, в силу чего удается
найти оптимальный компромисс для частотно-временного разрешения анализируемого сигнала.
Рис. 12 иллюстрирует различия в результатах применения методов FFT, STFT и CWT для анализа реакции усвоения ритмической
фотостимуляции в паттернах электроэнцефалограммы человека.
Амплитуда паттерна ЭЭГ до, во время и после фотостимуляции
(рис. 12, а), спектр мощности, построенный после быстрого преобразования Фурье этого паттерна (рис. 12, б), локальный вейвлетный
спектр |W(f, t0)|2паттерна (рис. 12, г) и нормированные интегральные распределения энергии вейвлетных спектров паттерна (сплошная линия рис. 12, д) и фотостимула (штрихпунктирная линия рис.
12, д); нормированные интегральные распределения энергии спектров, полученных оконным преобразованием Фурье для различных
значений ширины окна (w) и числа значений, которые перекрываются в каждом сегменте (n) (рис. 12 в, е), nL – общее число значений
в анализируемом паттерне.
24
г)
f (Гц)
20
0
8
10
Спектр
мощности
б)
14
8
0
5
10
f (Гц)
20
15
STFT, w = nL/5, n = w/2
1
0,5
8
10
12
t0 (с)
10
д)
2
0
г
12
16
FFT
4
в)
E(t0)/Emax(t0)
12
t (c)
14
16
E(t0)/Emax(t0)
–20
0
CWT
12,2
11,8
12
t0 (с)
14
16
14
16
CWT
1
0,5
0
е)
E(t0)/Emax(t0)
ЭЭГ (мкВ)
а)
8
10
12
t0 (с)
STFT, w = nL/9, n = w/6
1
0,5
0
8
10
12
t0 (с)
14
16
Рис. 12. Пример анализа реакции усвоения
ритмической фотостимуляции частотой 12 Гц
Рост амплитуды паттерна ЭЭГ, увеличение энергии локального
вейвлетного спектра |W(f, t0)|2и нормированного интегрального распределения энергии вейвлетного спектра в узком диапазоне частот
вблизи частоты фотостимуляции свидетельствуют о выраженной
реакции усвоения ритма внешней частоты. При этом быстрое преобразование Фурье позволяет обнаружить в спектре сигнала ЭЭГ
единственный максимум на частоте стимуляции (12 Гц), но не определить момент появления в сигнале этой частоты.
Как видно из рис. 12, нормированные интегральные распределения энергии спектра, полученного методом оконного преобразования Фурье, зависят от параметров w и n метода и, следовательно,
не позволяют единственным образом определить динамику удержания ритма заданной частоты и время запоминания ритма этой
частоты.
Непрерывное вейвлет преобразование, напротив, наглядно
показывает рост энергии локального вейвлетного спектра |W (f,
t0)|2(рис. 12, г) и его нормированного интегрального распределения
(сплошная линия на рис. 12, д) во время фотостимуляции и их спад
после выключения стимула.
25
3.8. Кросс-вейвлетные спектры
и оценивание вейвлет-когерентности двух сигналов
Для сравнения динамики двух сигналов y(t) и x(t) используется
оценивание кросс-вейвлетного спектра
Wxy = Wx (f, t0 )Wy* (f, t0 ) ,
который определяет локальные соотношения между двумя сигналами в определенные моменты времени на определенных частотах
[Maraun, Kurths, 2004].
На рис. 13 представлены два периодических сигнала x(t) и y(t),
заданные функциями:
 2 sin(2πw=
1t), w1 8, åñëè 0 < t < 2,2,
x(t) = 
πw2t), w2 2, åñëè 2á2 ≤ t < 4,6,
 0,5 sin(2=
y(t) = sin(2πw1t) + sin(2πw2t), w1 = 8, w2 = 2
и кросс-вейвлетный спектр этих сигналов, отражающий их ковариантность на частоте 8 Гц в первые 2,2 сек и последующую ковариантность на частоте 2 Гц в последние 2,4 сек.
Сигналы x(t) и y(t)
x(t), y(t)
4
2
0
–2
–4
0
0,5
1
1,5
2,5
2
3
3,5
4
4,5
t (c)
Кросс вейвлет спектр
Частота (Гц)
1
2
4
8
16
32
0
200
400
600
nL
800
1000
Рис. 13. Cигналы и их кросс-вейвлетные спектры
26
1200
Кроме того, для сравнения динамики двух сигналов y(t) и x(t)
можно оценить вейвлет-когерентность этих сигналов [Grinsted,
et al., 2004]:
WC =
S(Wxy )
S(Wxx )
S(Wyy )
.
Cимвол S обозначает предварительное сглаживание сигнала до
применения к нему вейвлет-преобразования с целью улучшения
частотно-временного разрешения и статистической значимости
[Maraun, Kurths, 2004].
Вейвлет – когерентность WC может принимать значения от 0 до 1
и описывает локальные корреляции между двумя сигналами; чем
ближе это значение к 1, тем более коррелированы сигналы.
Вопросы:
1. В чем состоит дискретное вейвлет- преобразование временного
ряда?
2. В чем состоит непрерывное вейвлет-преобразование временного ряда?
3. В чем различие между результатами непрерывного вейвлетпреобразования, быстрым преобразованием Фурье и оконным преобразованием Фурье?
4. Что такое локальный вейвлетный спектр энергии?
5. Что характеризует глобальный вейвлетный спектр энергии?
6. В чем различие между кросс-вейвлетным спектром и локальным вейвлетным спектром?
7.Как определить вейвлет-когерентность двух сигналов?
Лабораторная работа 1. Выполнить с помощью базисного вейвлета Добеши db4 дискретное вейвлет-преобразование временного
ряда и привести его графическую визуализацию
Лабораторная работа 2. Выполнить с помощью базисного вейвлета Морле непрерывное вейвлет-преобразование временного ряда и
привести его графическую визуализацию
Лабораторная работа 3. Применить вейвлет-анализ к исследованию биологических временных рядов (электроэнцефалограмма
человека в состоянии покоя и при гипервентиляции).
Лабораторная работа 4. Применить вейвлет-анализ к исследованию финансовых временных рядов (значений доходностей ценных бумаг).
27
Лабораторная работа 5. Применить вейвлет-анализ к исследованию объемов продаж на фондовых рынках.
Лабораторная работа 6. Применить вейвлет-анализ к исследованию динамики изменений в структуре паттернов ЭЭГ при действии внешнего периодического стимула.
Лабораторная работа 7. Применить вейвлет-анализ к подавлению шумовой составляющей сигнала.
28
4. ФРАКТАЛЬНЫЙ АНАЛИЗ ВРЕМЕННЫХ РЯДОВ
4.1. Понятие фрактальности временного ряда
Оценка степени фрактальности сигнала важна для исследования
флуктуаций физических, биологических и финансовых временных
рядов с целью получения информации о долговременных корреляциях и, соответственно, прогнозирования поведения этих рядов.
Фрактальность временного ряда связана с нерегулярностью с одной стороны и с некоторой повторяемостью в широком диапазоне
временных масштабов с другой стороны [Mandelbrot, 1982]. Для
функции, обладающей свойством фрактальности, справедливо следующее соотношение:
g(x0 + lx) − g(x0 ) ≈ l H (g(x0 + x) − g(x0 )),
где показатель степени H, называемый экспонентой Херста, характеризует нерегулярность функции g(x) в окрестности точки x0.
Чем меньше H, тем более нерегулярна функция [Павлов, Анищенко, 2007].
4.2. Оценивание степени фрактальности временного ряда
методом анализа флуктуаций относительно тренда
Для численной оценки показателя Херста применяется метод анализа флуктуаций относительно тренда (DFA – detrended fluctuation
analysis) [Peng, 1995].
Алгоритм вычисления показателя Херста состоит в следующем:
N
1) для исходного ряда значений {x(ti )}i =1 вычисляется интегриN
рованная последовательность {y(k)}k=1 , состоящая из накоплен
ных отклонений от среднего xk :
=
y(k)
k
∑ (xi − xˆk );
i =1
2) эта последовательность разбивается на некоторое число m = N/n
неперекрывающихся интервалов длины n;
3) в каждом из интервалов полученная интегрированная последовательность аппроксимируется прямой по методу наименьших
квадратов;
4) координата yn(k) прямой является средним значением интегрированной последовательности на k-ом интервале и обозначает
локальный тренд в пределах выбранного интервала;
29
5) среднеквадратичные отклонения (флуктуации) вычисленной интегрированной последовательности относительно локального тренда
для k-ого интервала вычисляются по формуле:
=
F (n)
1 N
∑ [y(k) − yn (k)]2 .
N k=1
6) Далее вычисления повторяются для других значений длины
интервала n от 5 до 100.
В силу того, что при увеличении длины интервала n значение
F(n), как правило, возрастает по степенному закону:
F(n) ~ nα,
скейлинговая экспонента α может быть вычислена как угловой коэффициент прямой, определяющей зависимость log F(n) от log n.
При α ≤ 1 полученное значение α совпадает со значением показателя Херста H.
Фрактальная размерность сигнала D связана с показателем Херста H следующим образом:
D = 2 – H.
К недостаткам метода следует отнести тот факт, что для надежного вычисления показателя Херста и, соответственно, фрактальной размерности, требуется большой временной ряд, а исследуемые
временные ряды, как правило, многократно меняют характер своего поведения при длительной записи. Поэтому вычисленная фрактальная размерность не будет связана с локальной динамикой изучаемого процесса, а будет характеризовать некоторую ее усредненную динамику. В этом смысле для описания динамики временного
ряда лучше подходит концепция мультифрактальности.
4.3. Оценивание степени мультифрактальности сигналов
методом поиска максимумов модулей вейвлет коэффициентов
Если для описания фрактальности монофрактального сигнала
достаточно одной величины (фрактальной размерности сигнала),
описывающей сохраняемость статистических характеристик исследуемого сигнала при изменении масштаба (частоты), то для того,
чтобы охарактеризовать мультифрактальный сигнал, необходимо
использовать целый спектр фрактальных размерностей.
Информация о возможной мультифрактальной структуре исследуемого сигнала и ее временной локализации t0 отражается в ассимптотическом поведении коэффициентов вейвлет-разложения сигнала W(f, t0) при малых значениях масштаба a, и, соответственно,
30
больших значениях частот f [Павлов, Анищенко, 2007]. Чем быстрее уменьшаются вейвлет коэффициенты при a → 0, тем более
регулярен сигнал в окрестности точки t0. Медленное уменьшение
вейвлет коэффициентов при a → 0 в окрестности точки t0 свидетельствует о наличии сингулярности (особенности в виде изрезанности
сигнала) в этой точке. Таким образом, скорость изменения модуля
коэффициентов вейвлет-преобразования сигнала позволяет определять наличие или отсутствие сингулярностей этого сигнала.
Степень сингулярности сигнала x(t) в точке t0 описывается экспонентой Гельдера, h(t0), наибольшей экспонентой, при которой
анализируемый сигнал может быть представлен суммой двух компонент: полинома Pn(t), описывающего регулярное поведение анализируемого сигнала и слагаемого, описывающего его нерегулярное
поведение [Павлов, Анищенко, 2007]:
x(t=
) Pn (t) + c t − t0
h(t0 )
.
В силу простой зависимости W (a,t0 ) ~ ah(t0 ) при a → 0 [Muzy et al.,
1993], экспонента Гельдера может быть вычислена следующим образом:
h(t0 ) ~
log10 W (a,t0 )
.
log10 a
При возрастании величины a влияние соседних сингулярностей может приводить к неточности вычисления, поэтому, как правило, экспонента Гельдера вычисляется на основании статистического описания
локальных сингулярностей с помощью частичных функций Z (q, a),
которые строятся по методу поиска максимумов модулей вейвлет коэффициентов (WTMM – wavelet transform modulus maxima) вдоль каждой
линии на масштабах, меньших заданного значения а [Muzy et al., 1993].
Детальное описание метода и его применение к экспериментальным
данным дано в работе [Arneodo et al., 1995].
Иллюстрация к поиску линий локальных максимумов проекции трехмерной вейвлетной поверхности |W(f, t0)|2сигнала, определяющей его локальный вейвлетный спектр, дана на рис. 14.
Алгоритм оценивания мультифрактальности сигнала состоит из
следующей последовательности процедур:
1) Применяется непрерывное вейвлет преобразование временного ряда, описывающего исследуемый сигнал x(t):
=
W
(a,t0 )
1
a
+∞
∫ x(t)ψ
−∞
*  t − t0


 dt.
 a 
31
f (Гц)
15
10
5
0
2
4
8
6
10
12
14
16
18
16
18
t0 (с)
Линии локальных максимумов
f (Гц)
15
10
5
0
2
4
6
8
10
12
14
t0 (с)
Рис. 14. Локальный вейвлетный спектр сигнала
и линии его локальных максимумов
2) Для каждого значения a находится множество L(a) линий локальных максимумов модулей вейвлет коэффициентов (линий, для
которых выполняется условие
∂ W (a,t0 )
= 0 ). Пример этого множе∂t
ства представлен на рис. 14.
3) Вдоль каждой линии для значений масштабов, меньших заданного значения a, вычисляются частичные функции Z(q, a) как
сумма q степеней максимумов модулей вейвлет коэффициентов
вдоль каждой линии на масштабах, меньших заданного значения а
(так называемые обобщенные статистические суммы):
Z (q, a) =
∑
l∈L(a)
(sup
a* ≤a
W (a* ,tl (a* ))
),
q
где tl(a*) определяет положение максимума, соответствующего линии l на этом масштабе.
32
τ(q)
4) В силу того, что при a → 0 частичная функция Z (q, a) ~ a
[Arneodo et al., 1995], скейлинговая экспонента τ(q) находится следующим образом:
τ(q) ~ log10 Z (q, a) log10 a.
4) Выбирая различные значения степенп q, можно получить линейную или нелинейную зависимость τ(q), что дает или постоянное
значение экспоненты Гельдера h(q) = const для монофрактальных
сигналов или большое число экспонент
h(q) =
dτ(q) dq ≠ const
для мультифрактальных сигналов, и получить, в этом случае, распределение экспонент Гельдера, называемое спектром сингулярности, которое вычисляется на основе скейлинговых экспонент:
D=
(h) qh(q) − τ(q).
4.4. Спектр сингулярности и его характеристики
Используя глобальные вейвлетные спектры и метод максимумов
модулей вейвлет преобразования, можно получить значения максимумов глобальной энергии тремора Emax и двух мультифрактальных параметров: 1) ширины спектра сингулярности
∆h = hmax – hmin
и 2) асимметрии спектра сингулярности
∆ = | ∆2 – ∆1 |,
где hmax= h (q = –5) и hmin= h (q = 5) – максимальное и минимальное
значения экспоненты Гельдера, соответствующие минимальной и
максимальной флуктуациям сигнала,
∆1 = hmax – h0 и ∆2 = h0 – hmin, h0 = h (q = 0).
Первый параметр Δh характеризует степень мультифрактальности анализируемого сигнала (чем больше величина Δh, тем больше
степень мультифрактальности).
Второй параметр ∆ определяет где, в области сильных сингулярностей (больших флуктуаций) (при q > 0), или в области слабых сингулярностей (малых флуктуаций) (при q < 0), в основном находятся
экспоненты Гельдера [Павлов, Анищенко, 2007].
Эти параметры определяют также степень коррелированности
интервалов в паттернах сигнала, так известно, что величины h < 0,5
соответствуют антикоррелированной динамике последовательных
значений сигнала (когда за большим значением сигнала следует ма33
лое значение), а величины h > 0,5 – его коррелированной динамике
(когда за большим значением сигнала следует также большое значение, а за малым – малое), поэтому в целом сигнал является более
гладким во втором случае, чем в первом [Xu et al., 2011].
На рис. 15 представлены зависимости скейлинговой экспоненты τ(q)
и экспоненты Гельдера h(q), полученные для тремора руки человека
двумя методами мультифрактального анализа (методом поиска максимумов модулей вейвлет коэффициентов (WTMM алгоритм) и методом мультифрактального анализа флуктуаций относительно тренда
(MF-DFA алгоритм)). Обе зависимости нелинейны, что свидетельствует
о большом числе экспонент Гельдера h(q) для исследуемого сигнала.
На этом рис. также изображены спектры сингулярности, D(h),
вычисленные с помощью двух мультифрактальных алгоритмов.
Спектр сингулярности D(h) также демонстрирует значительную
степень мультифрактальности анализируемого сигнала.
Для сравнения приведем на рис. 16 кривые для типично монофрактального сигнала, а именно движения броуновской частицы
0
1
τ(q)
h(q)
–2
0,5
–4
–6
–5
0
q
5
0
–5
0
q
5
D
1
0,5
0
0
0,5
h
1
Рис. 15. Примеры спектров скейлинговых экспонент, t(q)
и спектров сингулярности, D(h) для непроизвольных колебаний руки,
полученных с помощью WTMM алгоритма (сплошные кривые)
и MF-DFA алгоритма (штрихпунктирные кривые)
34
Монофрактал и зашумленный монофрактал
0
1
–1
D
τ(q)
–0,5
0,5
–1,5
–2
–5
0
q
5
0
q
5
0
0,5
0,55
0,6
h
0,65
0,7
h(q)
1
0,5
0
–5
Рис. 16. Спектры скейлинговых экспонент, t(q) и спектры
сингулярности, D(h) монофрактального сигнала (сплошная кривая)
и его зашумленной версии (штрихпунктирная кривая)
(сплошная кривая), и сигнала, полученного суперпозицией броуновского движения и нормально распределенной выборки случайных чисел (штрихпунктирная кривая).
Рассматриваемый монофрактальный сигнал характеризуется
единственной экспонентой Гельдера h(q) = 0,6, а спектр скейлинговых экспонент, t(q), представляет собой линейную функцию от степени q. Для зашумленной версии монофрактального сигнала мы наблюдаем небольшое отклонение экспоненты Гельдера от постоянной
величины (h(q) ≠ const). Спектр сингулярности D(h), представленный
на рис. 16, имеет в этом случае, в отличие от точки (единственного
значения для монофрактального сигнала) квадрaтичную форму,
которая характеризует флуктуации, соответствующие появлению
интервала значений экспонент Гельдера шириной Δh = 0,04. Однако малая ширина спектра сингулярности не позволяет думать о потери монофрактальности зашумленным сигналом. Следовательно,
несмотря на зашумление, сигнал, описывающий броуновское движение, как был, так и остается монофрактальным.
35
Вопросы:
1. Что характеризуют скейлинговые экспоненты?
2. В чем отличие метода DFA (анализа флуктуаций относительно тренда) от метода максимумов модулей вейвлет-преобразования
(WTMM)?
3. Как определить коррелированность динамики временного ряда?
4. Как связаны фрактальная размерность сигнала и показатель
Херста?
5. Чем отличаются монофракталы и мультифракталы?
6. Что характеризует спектр сингулярности?
7. Как определить степень мультифрактальности временного ряда?
Лабораторная работа 1. Определить наличие длительных корреляций в тестовых сигналах (белый шум, розовый шум и винеровский случайный процесс).
Лабораторная работа 2. Определить наличие долговременных
корреляций в финансовых временных рядах.
Лабораторная работа 3. Определить степень мультифрактальности заданного временного ряда методом максимумов модулей
вейвлет-преобразования (WTMM).
Лабораторная работа 4. Выполнить мультифрактальный анализ заданного паттерна электроэнцефалограммы человека.
Лабораторная работа 5. Провести мультифрактальный анализ корреляционных свойств последовательности кардиоинтервалов в электрокардиограмме человека до и после функциональной нагрузки.
36
5. РЕКУРРЕНТНЫЙ АНАЛИЗ ВРЕМЕННЫХ РЯДОВ
Другим перспективным методом исследования временных рядов
является метод рекуррентного анализа, позволяющий визуализировать определенные закономерности в рядах [Marwan et al., 2007].
5.1. Понятие рекуррентности временного ряда
В основе метода лежит построение и анализ рекуррентных диаграмм. Рекуррентная диаграмма является графическим представлением матрицы
(
)
Ri,j (m, ε) = Θ ε − yi − yj , i, j = 1, ..., N,
которая служит для отображения m-мерной фазовой траектории состояний y(t) на плоскость, в которой координатные оси являются
осями времени, а черная точка соответствует повторению состояния во времени i в некоторое другое время j.
В формуле, представленной выше, N – количество рассматриваемых состояний; ε – размер окрестности точки y, ||...|| – норма, m –
размерность пространства вложения, q – функция Хевисайда (принимающая значения 0 или 1). Таким образом
1, yi ≈ yj
Ri,j (m, ε) =
.
0, yi ≠ yj
При этом значения матрицы 1 или 0 отмечаются на рекуррентной диаграмме черной или белой точками и отражают рекуррентность или ее отсутствие с точностью до ε – ошибки. Рекуррентность
определяется как достаточная близость состояния yj состоянию
yi, то есть рекуррентными являются состояния yj, попадающие
в m-мерную окрестность с радиусом ε и центром в yi.
Фазовая траектория состояний y(t) получается из временного
ряда {x(t)} методом временных задержек [Takens, 1981].
5.2. Рекуррентная диаграмма
Рекуррентные диаграммы показывают структурные свойства
паттернов анализируемого сигнала и изменения в паттернах на протяжении записи. Отдельные изолированные точки на диаграмме
свидетельствуют о сильных флуктуациях, происходящих в анализируемом сигнале, то есть его стохастичности. Наличие диагональных линий, параллельных «линии идентичности», проходящей
под углом 45°, говорит о сходстве паттернов в различные времена.
37
Таким образом, появление длинных диагональных линий демонстрирует появление периодичности в исследуемом процессе, при
этом вертикальные расстояния между диагональными линиями соответствуют периодам колебаний. Если эти расстояния различны,
то колебания являются квазипериодичными. Возникновение в диаграмме вертикальных или горизонтальных линий говорит о появлении паттернов, которые не изменяются во времени или изменяются
очень медленно. Нерегулярное появление черных зон (кластеров),
соответствующих скоплениям вертикальных и горизонтальных
линий, а также белых зон, указывает на нерегулярность процесса
[Marwan et al., 2007].
Сравним рекуррентные диаграммы для периодического, зашумленного периодического, хаотического и стохастического временных рядов. Соответствующие сигналы и их рекуррентные диаграммы представлены на рис. 17 и рис. 18.
Рекуррентные диаграммы построены при величине размерности
вложения m = 3 и временной задержке d = 90 для периодического
и d = 15 для зашумленного периодического сигналов, d = 40 для стоПериод 1
Зашумленный периодический
2
2
1
0
0
–1
–2
–2
0
50
100
0
50
Стохастический
Хаотический
2
5
0
0
-2
0
10
20
Время (с)
30
40
–5
0
10
20
Время (с)
Рис. 17. Сигналы с различной структурой паттернов
38
100
Хаотический
2000
1500
1500
1000
1000
i
i
Период 1
2000
500
0
500
0
1000
0
2000
500
0
0
500
1000
1000
2000
Стохастический
1000
i
i
Зашумленный периодический
1000
0
500
0
0
500
1000
Рис. 18. Рекуррентные диаграммы для сигналов,
представленных на рис. 17
хастического процесса (броуновское движение) и d = 5 для хаотических колебаний системы Лоренца. Величина размера окрестности ε = 0,03σ для периодического сигнала, e = 2σ для зашумленного
периодического сигнала, e = 0,1σ для системы Лоренца и e =0,8σ для
броуновского движения, где σ – величина стандартного отклонения
анализируемых временных рядов.
Рекуррентная диаграмма периодического сигнала имеет длинные
непрерывающиеся диагональные линии с равными вертикальными
расстояниями между ними, соответствующими значению периода.
Рекуррентная диаграмма зашумленного периодического сигнала имеет равноотстоящие черные и белые полосы, соответствующие
зашумлению периода.
Стохастический процесс, описывающий движение броуновской частицы, имеет прерывающиеся нерегулярные белые и черные структуры.
Хаотическая динамика колебаний системы Лоренца приводит
к появлению коротких диагоналей и вертикальных структур, которые нерегулярны в отличие от периодического процесса. Други39
ми словами, хаотический процесс вызывает возникновение более
сложных квазипериодических рекуррентных структур с различными расстояниями между диагональными линиями, которые группируются в нерегулярные черные кластеры.
5.3. Численные показатели рекуррентных диаграмм
Используя метод рекуррентного анализа, можно определить численные показатели рекуррентных диаграмм [Marwan et al., 2007],
такие как:
1) показатель детерминизма (предсказуемости) исследуемого процесса
N
∑
=
DET
lP(ε, l)
l =lmin
N
−
∑ Ri,j (m,ε)
i,j
отношение числа рекуррентных точек, составляющих диагональ) {li , =
i 1, ...
ные структуры, к общему числу рекуррентных точек, где P(ε, l=
P(ε, l=
) {li , =
i 1, ..., Nl } – частотное распределение диагональных линий
длин l в рекуррентной диаграмме, N – число всех диагональных линий;
2) показатель дивергенции
DIV = 1/Lmax,
=
Lmax max
=
где
({li ,i 1,..., Nl }) – наибольшая диагональная линия
в рекуррентной диаграмме;
3) рекуррентное время
T=
j
{i, j : yi ,yj ∈ Ri ,yj −1 ∉ Ri } ,
необходимое для тогo, чтобы траектория вернулась в ε окрестность
точки, в которой она была ранее, определяемое как расстояние по
вертикали между началом и концом последовательной рекуррентной структуры в рекуррентной диаграмме, где Ri – рекуррентные
точки, которые принадлежат состоянию yi [Gao, 1999];
4) энтропия плотности рекуррентных времен
Tmax
EDRT = −
∑ P(T)ln P(T)
T =1
ln Tmax
,
где Tmax – наибольшее рекуррентное время в рекуррентной диаграмме, P(T) – функция плотности рекуррентных времен, полу40
ченная после нормирования гистограммы H(T) рекуррентных
времен:
P(T) = T
H(T)
max
.
∑ H(T)
T =1
Показатель детерминизма DET, или предсказуемости поведения
процесса, максимален (= 1) для периодического процесса и минимален (= 0) для стохастического, так как процессы со стохастическим
поведением порождают очень короткие диагонали, либо вообще не
имеют их, в то время как детерминистские процессы дают длинные
диагонали [Marwan et al., 2007]. Показатель дивергенции, наоборот, минимален (= 0) для периодического процесса и характеризует
сходимость сегментов траектории (то есть чем дальше процесс от периодического, тем короче диагональные линии и тем больше значение показателя DIV) [Marwan et al., 2007].
Гистограмма рекуррентных времен H(T) имеет единственный
максимум для чисто периодического сигнала и множество максимумов для сильно зашумленного неоднородного сигнала, а функция плотности рекуррентных времен P(T) принимает значения от 1
для периодического сигнала до 1/Tmax для неоднородного сигнала.
Энтропия плотности рекуррентных времен EDRT отражает сложность детерминистской составляющей сигнала и характеризует неопределенность его периода. Эта величина может быть равной 0 (что
означает отсутствие неопределенности в значении периода для строго
периодического сигнала) или равной 1 (максимальная неопределенность для стохастического сигнала) [Little et al., 2007]. Для реальных
сигналов EDRT принимает промежуточные от 0 до 1 значения.
Таким образом, метод рекуррентного анализа, позволяющий не
только визуализировать определенные закономерности в паттернах
сложных нестационарных колебаний, но и выявлять количественные
параметры эволюции сигнала во времени, является перспективным
методом исследования нелинейной динамики сложных сигналов.
5.4. Гистограмма плотности рекуррентных времен
Построим для хаотических колебаний системы Лоренца гистограмму функции плотности рекуррентных времен, необходимых
для возврата траектории в ε окрестность точки, в которой она была
ранее (рис. 19). Штрихпунктирной линией обозначена гистограмма,
вычисленная для суррогатных данных (рандомизированных коле41
Сигнал зашумленной системы Лоренца
Плотность рекуррентных времен
1
0,9
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
0
2
8
10
Рис. 19. Гистограмма функции плотности рекуррентных времен
для колебаний системы Лоренца (сплошная кривая)
и для суррогатного сигнала, полученного рандомизацией исходных
колебаний (штрихпунктирная кривая). Размерность вложения m = 3,
временная задержка d = 5 и окрестность e = 0,1 σ
баний системы Лоренца). Максимум гистограммы рекуррентных
периодов для колебаний системы Лоренца равен 5.75 ± 0.3 (с).
5.5. Локализация неустойчивых периодических орбит
Паттерн, соответствующий периодическим колебаниям (периодическим орбитам), отражается в рекуррентной диаграмме непрерывающимися, равноотстоящими друг от друга диагональными
линиями. Вертикальные расстояния между диагональными линиями соответствуют периодам колебаний. Хаотические паттерны
приводят к появлению коротких диагоналей, вертикальные расстояниями между которыми становятся нерегулярными.
Когда траектория колебательного процесса приближается к неустойчивой периодической орбите, она остается в его окрестности
на некотором временном интервале, длина которого зависит от
того, насколько неустойчива орбита [Dhamala, et al., 2000]. Поэтому неустойчивые периодические орбиты могут быть обнаружены
путем поиска внутри рекуррентной диаграммы окон, в которых
паттерны соответствуют периодическому движению. Если расстояния между диагональными линиями изменяются от одного
42
выбранного окна к другому, то существуют орбиты с различными
периодами. Период неустойчивой периодической орбиты можно
оценить по вертикальному расстоянию между рекуррентными
точками в окне, деленному на частоту дискретизации временного
ряда [Dhamala, et al., 2000].
Алгоритм нахождения неустойчивой периодической орбиты состоит из следующих процедур.
1. Поиск рекуррентных времен с точностью до ε – ошибки, то
есть в окрестности с радиусом ε.
2. Значения рекуррентных периодов определяются как рекуррентные времена, деленные на частоту дискретизации. На основании этих значений строится гистограмма рекуррентных периодов.
Значения периодов неустойчивых периодических орбит находятся
как максимумы этой гистограммы.
3. Чтобы исключить влияние шума, найденные неустойчивые
периодические орбиты проверяются на статистическую значимость. Для этого описанные выше процедуры повторяются для не
менее 30 суррогатных сигналов, полученных рандомизацией исходных данных [Theiler, et al., 1992].
Тремор, период 6
2
15
1,5
10
1
5
0,5
x(t + d)
x(t + d)
Система Лоренца, период 1
20
0
0
–5
–0,5
–10
–1
–15
–1,5
–20
–20
–10
0
10
20
–2
–2
–1
0
x(t)
1
2
Рис. 20. Неустойчивая орбита периода 1 для колебаний системы Лоренца
и орбита периода 6 для непроизвольных колебаний руки человека
43
Статистическая мера наличия статистически значимых неустойчивых периодических орбит в исследуемом сигнале определяется
отношением
_
k =( A − A) / σ,
_
где A – значение максимума гистограммы, A – значение A для
суррoгатных сигналов и σ – стандартное отклонение. Величина k
характеризует существование статистически значимых неустойчивых периодических орбит в исследуемом сигнале по сравнению с его
рандомизированным вариантом. Значение k > 2 означает обнаружение неустойчивых периодических орбит с более чем 95% уровнем
значимости [Bevington, et al., 2003].
Орбита периода 1, найденная для хаотических колебаний системы Лоренца, и орбита периода 6 для непроизвольных колебаний
руки человека показаны на рис. 20.
5.6. Кросс-рекуррентный анализ
и построение совместной рекуррентной диаграммы
Метод кросс-рекуррентного анализа позволяет визуализировать
определенные закономерности в двух сигналах (даже от различных
физических источников) [Marwan, et al., 2007].
Совместная рекуррентная диаграмма (СРД) является графиче1, yi ≈ yj , xi ≈ xj ,
ским представлением матрицы Ri,j (ε) =
, в которой
 0, yi ≠ yj xi ≠ xj
значения 1 или 0 соответствуют черной или белой точкам, при этом
черная точка соответствует рекуррентности, а белая – ее отсутствию.
Совместная рекуррентность, с точностью до ε – ошибки, определяется как возврат состояния yj к состоянию yi и одновременный возврат
состояния xj к состоянию xi. Величина ε обычно выбирается равной 1%
от величины стандартного отклонения анализируемого сигнала.
К численным показателям совместных рекуррентных диаграмм относятся такие показатели как: L – средняя длина диагональных линий
в СРД и τ – рекуррентное время, необходимое для тогo, чтобы значение
сигнала вернулось вε окрестность точки, в которой оно было ранее.
На рис. 21 представлены примеры совместных рекуррентных диаграмм паттернов ЭЭГ и ритмических световых сигналов с частотой
12 Гц для пациентов с мерцательной аритмией двух типов (пароксизмального и постоянного, то есть различающихся по длительности
нарушений сердечного ритма). Эти диаграммы построены при величине временной задержки d = 5 и размерности вложения m = 3,
44
Паттерны ЭЭГ и световые сигналы,12 Гц
в)
40
а)
40
20
20
0
0
–20
–20
600
800
1000
1200
800
600
1000
1200
Совместные рекуррентные диаграммы
б)
г)
i
1000
i
1000
500
0
500
0
500
1000
j
0
0
500
1000
j
Рис. 21. Примеры паттернов ЭЭГ во время фотостимуляции пациентов
с мерцательной аритмией: а – пароксизмальной формы;
в – постоянной формы; б, г – совместные рекуррентные диаграммы
этих паттернов и световых сигналов
величина размера окрестности ε равна 1% от величины стандартного
отклонения анализируемых временных рядов.
Левая рекуррентная диаграмма содержит квазипериодические
рекуррентные структуры с различными расстояниями между диагональными линиями, которые группируются в нерегулярные черные кластеры, что свидетельствует о квазипериодическом возникновении одновременных рекуррентностей в паттерне ЭЭГ у пациента с мерцательной аритмией пароксизмальной формы и световом
сигнале заданной частоты.
В отличие от этого, правая диаграмма содержит только редкие и
короткие диагональные линии, то есть почти не имеет совместных
рекуррентностей в анализируемом паттерне пациента с фибрилляцией предсердий в постоянной форме и данном световом сигнале.
45
Средняя длина диагональных линий для левой диаграммы L = 7,8,
для правой – L = 2,1, рекуррентное τ = 0,66 и τ = 0,28 соответственно.
Таким образом, рекуррентный анализ позволяет обнаружить
достоверные различия в количественных параметрах реакции усвоения ритма фотостимуляции у разных групп пациентов с мерцательной аритмией, различающихся по времени существования нарушения сердечного ритма.
Вопросы:
1. Что такое рекуррентная диаграмма и чем отличаются рекуррентные диаграммы для заданного временного ряда и его зашумленного варианта?
2. Что определяет показатель детерминизма и как его определить?
3. Что определяет показатель дивергенции и как его определить?
4. Что такое рекуррентное время и как его определить?
5. Как определяется энтропия плотности рекуррентных времен?
6. Как построить гистограмму плотности рекуррентных времен
и что она характеризует?
Лабораторная работа 1. Построить рекуррентные диаграммы
для заданного временного ряда и его зашумленного варианта.
Лабораторная работа 2. Определить показатели рекуррентных
диаграмм (показатель детерминизма, показатель дивергенции, рекуррентное время) для заданного временного ряда и его зашумленного варианта.
Лабораторная работа 3. Построить гистограмму плотности рекуррентных времен и вычислить энтропию плотности рекуррентных
времен для заданного временного ряда и его зашумленного варианта.
46
6. ЗАДАЧА ОПТИМИЗАЦИИ ПОРТФЕЛЯ ЦЕННЫХ БУМАГ
6.1. Характеристика портфельных инвестиций
Основная идея стратегии инвестирования состоит в снижении
уровня риска путем выбора наилучшей комбинации ценных бумаг
(создания оптимального портфеля инвестиций). Оптимальность
портфеля определяется обеспечением заданного уровня дохода и
минимизацией риска.
Выделяют три типа портфелей (Касимов, 1998):
1) портфели роста: а) портфели агрессивного роста (формируются
из ценных бумаг быстрорастущих компаний, вложения в которые
являются рискованными, но ориентированными на максимальный
прирост капитала), б) портфели консервативного роста (формируются с целью сохранения инвестиционного капитала и незначительного его увеличения без существенного риска),
2) портфели дохода (формируются для обеспечения высокого
уровня текущего дохода в виде процентных выплат),
3) смешанные портфели (роста и дохода), которые комбинируют
свойства обоих портфелей.
Если в момент времени t ценная бумага имеет доходность ri, то
доходность портфеля, состоящего из ценных бумаг с долями (x1,
x2, …, xn), где
n
n
i
i
∑ xi = 1 определяется следующим образом R = ∑ xiri .
6.2. Оптимизация портфелей ценных бумаг
Cтандартное отклонение доходности портфеля представляет собой волатильность доходности портфеля, включающую в себя изменения доходности как в положительную, так и в отрицательную
стороны.
Задача нахождения оптимального портфеля состоит в нахождении минимума максимальных потерь, то есть минимума среднеквадратичных отклонений от средней доходности [Markowitz, 1959,
Касимов, 1998.].
С точки зрения математики, выбор оптимального портфеля может быть сделан на основе решения задачи квадратичного программирования, в которой осуществляется поиск максимума доходности портфеля
n
R = ∑ xiri
i
47
и минимума дисперсии доходности портфеля
=
V
n n
∑∑ xi xj σij ,
i
j
где sij – ‘элементы матрицы ковариации.
Эта задача является задачей нахождения условного экстремума
функции нескольких переменных и решается с помощью построения функции Лагранжа
m
L(x,=
λ) V (x) − ∑ λ i gi (x)
i
от m + n переменных и определения знака матрицы вторых частных
производных от функции Лагранжа, то есть матрицы Гессе с эле∂2 L
ментами Hij =
, где gi(x) = 0 – условия ограничений на доход∂xi ∂λ j
ность ценных бумаг.
Вопросы:
1. Что такое доходность портфеля и волотильность доходности
портфеля?
2. Как оптимизировать портфель ценных бумаг?
Лабораторная работа 1. Методом квадратичного программирования решить задачу оптимизации заданного портфеля ценных
бумаг.
48
Список литературы
1. Лоскутов А. Ю., Михайлов А. С. Теория сложных систем. Ижевск,
2007. 620 c.
2. Ivanov P. C., Amaral L. A., Goldberger A. L., et al. Multifractality
in human heartbeat dynamics. Nature. 1999. P. 461–465.
3. Безручко Б. П., Смирнов Д. А. Математическое моделирование и хаотические временные ряды. Гос. Саратов: УНЦ Колледж, 2005. 320 с.
4. Kantz H., Schreiber T. Nonlinear Time Series Analysis. Cambridge: Cambridge Univ. Press, 2002.
5. Takens F. Detecting strange attractors in turbulence. Dynamical Systems and Turbulence, Lecture Notes in Mathematics (D. Rand,
L. S. Young, eds.). Berlin: Springer-Verlag, 1981. Vol. 898. P. 366–381.
6. Fraser A. M., Swinney H. L. Independent coordinates for strange
attractors from mutual information. Phys. Rev. 1986. P. 1134–1140.
7. Kennel M. B., Brown R., Abarbanel H. D. Determining embedding
dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A. 1992. P. 3403–3411.
8. Астафьева Н. М. Вейвлет-анализ: основы теории и примеры
применения // Успехи физических наук. 1998. Т. 166. С. 1145–1170.
9. Дремин И. М., Иванов О. В., Нечитайло В. А. Вейвлеты и их использование // Успехи физических наук. 2001. Т. 171. № 5. С. 465 –501.
10. Ширяев В. И. Финансовые рынки и нейронные сети: учеб. пособие. M.: ЛКИ, 2007. 234 с.
11. Percival D. B., Walden A. T. Wavelet methods for time series
analysis. Cambridge: Cambridge University Press, 2000.
12. Малла С. Вейвлеты в обработке сигналов. М.: Мир, 2005. 671 с.
13. Добеши И. Десять лекций по вейвлетам: Регулярная и стохастическая динамика. Ижевск, 2001. 464 c.
14. Смоленцев Н. К. Основы теории вейвлетов. Вейвлеты в Matlab.
М.: ДМК Пресс, 2005. 301 с.
15. Блаттер К. Вейвлет-анализ. Основы теории. М.: Техносфера,
2006. 280 с.
16. Мусалимов В. М., Дик О. Е., Тюрин А. Е. Параметры действия
энергетического спектра вейвлет-преобразований // Известия вузов. Приборостроение. 2009. Т. 52. № 5. С. 10–15.
17. Grossmann A., Morlet J. Decomposition of Hardy functions into
square integrable wavelets of constant shape // S.I.A.M. Journal of
Mathematical Analysis. 1984. Vol. 15. P. 723–736.
18. Hu J., Gao J. B., Wang X. S. Multifractal analysis of sunspot
time series: the effects of the 11-year cycle and Fourier truncation.
J Stat Mech. 2009. P. 1–20.
49
19. Allen J. B. Short Time Spectral Analysis, Synthesis and Modification by Discrete Fourier Transform. IEEE Transactions on Acoustics, Speech and Signal Processing. 1977. P. 235–238.
20. Maraun D., Kurths J. Cross wavelet analysis: significance testing and pitfalls. Nonlinear Processes in Geophysics. 2004. P. 505–514.
21. Grinsted A., Moor J. C., Jevrejeva S. Application of the cross wavelet transform and wavelet coherence to geophysical timeseries. Nonlinear processes in geophysics, 2004, 11: 561–566.
22. Mandelbrot B. B. The fractal geometry of nature. San Francisco:
W. H. Freeman, 1982.
23. Павлов А. Н., Анищенко В. С. Мультифрактальный анализ
сложных сигналов // Успехи физических наук. 2007. Т. 177. С. 859.
24. Peng C. K., Havlin S., Stanley H. E., Goldberger A. L. Quantification of scaling exponents and crossover phenomena in nonstationary
heartbeat time series // Chaos: An Interdisciplinary Journal of Nonlinear Science. 1995. 5(1). 82–87.
25. Muzy J. F., Bacry E., Arneodo A. Multifractal formalism for fractal signals: the structure-function approach versus the wavelet-transform modulus-maxima method. Phys. Rev. E, 1993. 47: 875–884.
26. Arneodo A., BacryE., Muzy J. F. The thermodynamics of fractals
revisited with wavelets. Physica A. 1995. P. 232–275.
27. Xu Y., Ma Q. D. Y., Schmitt, D. T., et. al. Effects of coarse-graining on the scaling behavior of long-range correlated and anti-correlated signals. Physica A. 2011. P. 4057–4072.
28. Marwan N., Romano M. C., Thiel M., et. al. Recurrence plots for
the analysis of complex systems. Physics Reports. 2007. P. 237–329.
29. Gao J. B. Recurrence time statistics for chaotic systems and
their applications. Phys. Rev. Lett. 1999. P. 3178–3181.
30. Little M. A., McSharry P. E., Roberts S. J., et al. Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection. BioMedical Engineering OnLine, 2007. Vol. 6.
31. Dhamala M., Lai Y.-C., Kostelich E. J. Detecting unstable periodic
orbits from transient chaotic time series. Phys. Rev. E. 2000. P. 6485–6489.
32. Theiler J., Eubank S., Longtin A., et al. Testing for nonlinearity
in time series: the method of surrogate data. Physica D. 1992. P. 77–94.
33. Bevington J. R., Robinson D. K. Data reduction and error analysis for the physical sciences. Third edition. McGraw-Hill Higher Education, 2003. 338 p.
34. Касимов Ю. Ф. Основы теории оптимального портфеля ценных бумаг. М.: Филин. 144 c.
35. Markowitz H. M. Portfolio Selection: Efficient Diversification of
Investment. New York: Wiley, 1959.
50
СОДЕРЖАНИЕ
1. Понятие временного ряда, динамической системы и аттрактора
динамической системы.............................................................
1.1. Свойства временных рядов: стационарность,
нестационарность, стохастичность, детерминизм................... 1.2. Понятие аттрактора динамической системы.................... 1.3. Реконструкция аттрактора из исходного временного ряда....
2. Вычисление корреляционного интеграла и корреляционной
размерности восстановленного аттрактора...................................
2.1. Алгоритм вычисления корреляционной размерности
аттрактора........................................................................ 3. Вейвлет-анализ временных рядов...........................................
3.1. Дискретное вейвлет преобразование временного ряда........ 3.2. Интегральная оценка спектральной плотности энергии
временного ряда................................................................. 3.3. Непрерывное вейвлет преобразование временного ряда
и оценивание энергии сигнала.............................................. 3.4. Базисный вейвлет Морле............................................... 3.5. Глобальный вейвлетный спектр..................................... 3.6. Удаление шумовой составляющей сигнала....................... 3.7. Сравнение вейвлетного преобразования сигнала
с быстрым преобразованием Фурье и оконным
преобразованием Фурье...................................................... 3.8. Кросс-вейвлетные спектры и оценивание
вейвлет-когерентности двух сигналов................................... 3
3
4
4
11
11
15
15
17
19
19
19
22
23
26
4. Фрактальный анализ временных рядов...................................
4.1. Понятие фрактальности временного ряда........................ 4.2. Оценивание степени фрактальности временного ряда
методом анализа флуктуаций относительно тренда................. 4.3. Оценивание степени мультифрактальности сигналов
методом поиска максимумов модулей вейвлет коэффициентов.....
4.4. Спектр сингулярности и его характеристики ................... 29
29
5. Рекуррентный анализ временных рядов..................................
5.1. Понятие рекуррентности временного ряда....................... 5.2. Рекуррентная диаграмма.............................................. 5.3. Численные показатели рекуррентных диаграмм .............. 5.4. Гистограмма плотности рекуррентных времен................. 5.5. Локализация неустойчивых периодических орбит............ 5.6. Кросс-рекуррентный анализ и построение совместной
рекуррентной диаграммы.................................................... 37
37
37
40
41
42
6. Задача оптимизации портфеля ценных бумаг...........................
6.1.Характеристика портфельных инвестиций ...................... 6.2. Оптимизация портфелей ценных бумаг........................... 47
47
47
Список литературы..................................................................
49
29
30
33
44
51
Учебное издание
Дик Ольга Евгеньевна
НЕЛИНЕЙНЫЙ АНАЛИЗ
ВРЕМЕННЫХ РЯДОВ
Учебное пособие
Публикуется в авторской редакции
Компьютерная верстка В. Н. Костиной
Сдано в набор 27.12.17. Подписано к печати 12.02.18. Формат 60 × 84 1/16.
Усл. печ. л. 3,0. Уч.-изд. л. 3,3.
Тираж 50 экз. Заказ № 57.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Документ
Категория
Без категории
Просмотров
3
Размер файла
6 582 Кб
Теги
dik
1/--страниц
Пожаловаться на содержимое документа