CC BY
143
УДК 519.2:004.9 Л. О. КИРИЧЕНКО
ВЕЙВЛЕТ-ОЦЕНИВАНИЕ ПОКАЗАТЕЛЯ ХЕРСТА ПО ВРЕМЕННЫМ РЯДАМ СО ЗНАЧИТЕЛЬНЫМИ ЦИКЛИЧЕСКИМИ КОМПОНЕНТАМИ С ПРИМЕНЕНИЕМ ПАКЕТНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ
Предлагается метод вейвлет-оценивания показателя Херста, основанный на использовании пакетного преобразования. Применение разработанного метода позволяет определять степень самоподобия для рядов небольшой длины, обладающих существенными трен-довыми и циклическими компонентами.
1. Введение
В настоящее время стало общепризнанным, что многие информационные, биологические, физические, технологические процессы имеют сложную фрактальную структуру. Фрактальный анализ используется для моделирования, анализа и контроля сложных систем в различных областях науки и техники. В частности, в геологии - для прогнозирования сейсмической активности и цунами, определения возраста геологических пород; в биологии и медицине - для диагностики заболеваний и физиологического состояния по записям ЭКГ и ЭЭГ, при исследовании мутаций и изменений на генетическом уровне; в экономике - для прогнозирования кризисных ситуаций и оценивания риска по финансовым рядам; в физике -для исследования турбулентности и термодинамических процессов [1-4].
Оценка степени самоподобия (показателя Херста) по экспериментальным данным играет важнейшую роль в изучении процессов, имеющих фрактальные свойства. Степень самоподобия стохастического процесса, который определяется показателем Херста, является показателем сложности динамики и корреляционной структуры временного ряда.
Для оценки показателя Херста по временному ряду существует ряд методов [5-6]. Однако большинство методов применимы только к стационарным рядам, в то время как множество природных, технических и информационных процессов являются нестационарными. Типичным нарушением стационарных свойств, которые встречаются на практике, является наличие трендовой и циклической составляющей ряда. Эффективным инструментом исследования структуры нестационарных временных рядов является аппарат вейвлет-анализа, который позволяет выделять структурные компоненты ряда. В частности, методы, основанные на дискретном вейвлет-преобразовании, позволяют оценивать степень самоподобия для нестационарных рядов [7-9].
Однако стандартный метод вейвлет-оценивания не позволяет получить адекватные оценки показателя Херста при наличии значительных циклических компонент временного ряда [10, 11]. Целью представленного исследования является разработка метода оценивания показателя Херста для рядов с существенными трендовыми и циклическими компонентами, основанного на применении пакетного вейвлет-преобразования.
2. Декомпозиция временных рядов с помощью дискретного вейвлет-
преобразования [7-9,12]
Вейвлет-преобразование временного ряда - это его представление в виде обобщенного ряда или интеграла по системе базисных функций
полученных из материнской вейвлет-функции у(1) за счет операций сдвига во времени ь и изменения временного масштаба а . В дискретном вейвлет-преобразовании (ДВП) параметры масштаба и сдвига задаются обычно в виде степенных функций: а = 2-т , ь = к0 -2-т , т, к0 е Ъ . Дискретные вейвлеты используются, как правило, в паре со связанными с ними дискретными скейлинг-функциями. Скейлинг-функции имеют с вейвлетами общую об-4
ласть задания и определенное соотношение между значениями. При заданных материнском вейвлете уОО и соответствующей скейлинг-функции ф^) аппроксимирующие коэффициенты аргЦк) и детализирующие коэффициенты det(j,k) ДВП для процесса Х^) определяются следующим образом:
да да
аргЦк) = | Х(0фdet(j,k) = | X(t)уj,k(t)dt,
—да —да
где j - параметр масштаба; к - параметр сдвига (j, к е 2) и базисные функции равны:
ф ^а)=2—j/2 ф( 2—jt—к),
у j,k (t) = 2—j/2 у(—jt — к).
Декомпозиция временного ряда, выполняемая с помощью ДВП, заключается в разбиении исследуемого ряда на две составляющие: аппроксимирующую и детализирующую, с последующим аналогичным дроблением аппроксимирующей до заданного уровня разложения. Временной ряд X(t) представляется в виде суммы аппроксимирующей компоненты арргохм и детализирующих компонент detailj (t):
N
Ха)=арргохм (t)+^ detailj (t) =
j=l
Na N Nj
= Е apr(N,k) фJ)k(t) + Е Е det(j, к)уj,k(t), к=1 j=1 к=1
где N - выбранный максимальный уровень разложения; N - количество детализирующих
коэффициентов на уровне j; Na - количество аппроксимирующих коэффициентов на уровне
N .
Выбор типа вейвлет-функции и количества уровней разложения является важным вопросом при выполнении ДВП. Обычно вейвлет-функция подбирается в зависимости от временных и частотных характеристик каждого анализируемого сигнала. Максимальный уровень разложения зависит от того, какие частотные диапазоны необходимо исследовать.
В обычном алгоритме ДВП при переходе с масштабного уровня т на уровень т+1 функция аппроксимирующих коэффициентов apr(m,k) разделяется на низкочастотную apr(m + 1,к) и высокочастотную det(m + 1,к) части спектрального диапазона. При дальнейшем увеличении масштабных уровней разложению подвергаются только низкочастотные аппроксимирующие функции, детализирующие коэффициенты далее не анализируются (рис. 1 слева)).
Рис. 1. Дерево обычного ДВП (слева) и пакетного ДВП (справа)
Идея пакетного вейвлет-анализа состоит в том, чтобы раскладывать детализирующие функции, используя тот же самый метод разложения. В пакетном алгоритме ДВП операция последовательного частотного расщепления применяется как для низкочастотных аппроксимирующих, так и для высокочастотных детализирующих коэффициентов. Таким образом пакетное ДВП способствует лучшей частотной локализации. В результате возникает дерево разложения, пример которого показан на рис. 1 справа. При разложении данных некоторого узла ниже слева находятся аппроксимирующие коэффициенты, а справа -детализирующие.
Каждый узел дерева пакетного ДВП содержит набор вейвлет-коэффициентов, соответствующих определенному частотному диапазону. Временной ряд может быть однозначно восстановлен по пакетным вейвлет-коэффициентам, находящимся в терминальных узлах дерева.
3. Оценивание самоподобия стохастических процессов с помощью ДВП [7-10]
Стохастический процесс Х(г) является самоподобным с параметром самоподобия н, если процесс а-нХ (аг) описывается теми же конечномерными законами распределений,
что и Х(г): Law{X(t)} = Law{a-нХ(аг)}, va > 0,1 > 0 .
Параметр Н, 0 < Н < 1, называемый показателем Херста, представляет собой степень самоподобия и меру долгосрочной зависимости стохастического процесса.
Метод оценивания степени самоподобия базируется на свойствах детализирующих вейвлет-коэффициентов, полученных в результате декомпозиции процесса. Если случайный процесс Х(г) является самоподобным, то детализирующие коэффициенты на каждом уровне разложения также обладают свойством самоподобия:
Law{det(j,k)} = Law{2j(H+1/2) ай(0,к)}, (1)
где det(j,k) - к-й детализирующий коэффициент уровня разложения j, к = 1,2,...,^; Nj -
количество вейвлет-коэффициентов на уровне разложения j; Н - показатель Херста. Это свойство вытекает из свойств вейвлет-базиса, полученного путем масштабирования материнских вейвлетов.
Если для самоподобного процесса Х(г) существуют моменты 2-го порядка, то для вейвлет-коэффициентов, полученных в результате декомпозиции данного процесса, выполняется следующее равенство:
М^ак)|2 = М^е1(0,к)|2 2"](2Н+1). (2)
Величина вейвлет-энергии М^егЦк)2 на уровне j вычисляется следующим образом:
1 N
— Е лл'аю. (3)
Ч к=1
Метод вейвлет-оценивания степени самоподобия Н базируется на свойствах детализирующих вейвлет-коэффициентов, полученных на разных уровнях декомпозиции исследуемого временного ряда. Основу данного метода составляет положение, следующее из (1)-(2), о том, что изменение значений вейвлет-энергии на масштабном уровне j подчиняется скейлинговому отношению:
Е х 2(2Н+1)-'. (4)
График зависимости 1оя2(Ер от j представляет собой прямую линию, аппроксимированную методом наименьших квадратов. Значение параметра н можно найти, оценив тангенс угла наклона прямой. На рис. 2 представлен типичный спектр вейвлет-энергии
самоподобного процесса (вверху) и зависимость 1°Я2(Ер от номера уровня разложения j (внизу).
0.4
0.2
0 0
Рис. 2. Спектр вейвлет-энергии самоподобного процесса (вверху) и зависимость 1^2 (Ер (внизу) Одной из важных характеристик вейвлет-функции у^) является наличие нулевых мо-
ментов высокого порядка: | ^у(^ = 0, к = 0 ,р 1. Число нулевых моментов р характе-
—да
ризует способность вейвлет-преобразования анализировать мелкомасштабные флуктуации и особенности высокого порядка, игнорируя регулярные полиномиальные составляющие степени т < р .
4. Моделирование реализаций фрактального процесса
Одной из наиболее известных и простых моделей стохастической динамики, обладающих фрактальными свойствами, является фрактальное броуновское движение (ФБД), которое широко применяется в физике, химии, биологии, экономике и теории сетевого трафика [1-4].
Гауссовский процесс Х^) называется фрактальным броуновским движением с параметром Н, 0 < Н < 1, если приращения случайного процесса ДХ(х) = Х^ + т) — Х^) имеют распределение вида:
Р(ДХ < x) = -
1
J ExP
2CTq Т
dz
л/2лсто тН
где сто - коэффициент диффузии.
ФБД с параметром Н = 0,5 совпадает с классическим броуновским движением. Приращения ФБД называются фрактальным гауссовским шумом (ФГШ), дисперсия которого
подчиняется соотношению Б[Х^ + т) — Х^)] = ст02т2Н .
Существует несколько методов построения ФБД для случая дискретного времени. Эти модели имеют определенные недостатки, одним из которых является занижение (завышение) степени самоподобия процесса при малых или больших теоретических значениях параметра Херста и небольшой длине модельной реализации [13].
Одним из методов, направленных на преодоление указанных недостатков, является построение ФБД с помощью биортогональных вейвлетов [14-15]. В этом случае построение реализации ФБД осуществляется с помощью дискретного вейвлет-преобразования, где детализирующие вейвлет-коэффициенты на каждом уровне являются независимыми нормально распределенными случайными величинами, а аппроксимирующие вейвлет-коэффи-циенты получены с помощью фрактального процесса авторегрессии и скользящего среднего БАММА:
да
Вн« = Е фна—к4Н) +
к=—да
да да
+Е Е 2—->н тн(2^—к)Е^ — Ъ0'
]=0к=—да
0.8
E
0.6
2
3
4
5
6
7
8
2
3
4
5
6
7
8
+<ю
где Тн - биорогональная базовая вейвлет-функция; Фн - соответствующая скей-линг-функция; skH) - стационарный гауссовский процесс FARIMA с параметром дробного дифференцирования d = H - 0,5; е j,k - независимые стандартные нормальные случайные
величины; bo - константа, такая что Вн(0) = 0 .
На рис.3 приведены реализация ФБД (вверху) и соответствующий ряд ФГШ (внизу), полученные описанным выше методом. Длина реализации равна 1024 значениям, теоретический показатель Херста Н=0.8.
Fractional Brownian motion
0 100 200 300 400 500 600 700 800 900 Fractional Gaussian noise
0 100 200 300 400 500 600 700 800 900 1000
Рис.3. Реализация ФБД (вверху) и ФГШ ( внизу)
5. Оценивание самоподобия временных рядов с циклическими составляющими
В работе использовались различные модельные сигналы - временные реализации, представленные суммой ФГШ с заданным показателем Херста и периодической составляющей. Суммарный сигнал X(t) можно представить как
X(t) = P(t) + fgn(t),
где P(t)- циклический тренд; fgn(t)- фрактальный шум. Соотношение циклического тренда и модельного фрактального шума можно определить как
Ratio =
S '
^noise
где Strend - среднее квадратическое отклонение тренда; Snoise- среднее квадратическое фрактального процесса.
Для оценивания показателя Н были выбраны вейвлет-функции семейства Добеши (db2-db20), которые обладают необходимым числом нулевых моментов (соответственно 2-20). При небольших значениях отношения тренда к ФГШ (Ratio < 2) метод вейвлет-оценивания (4) позволяет получить адекватные оценки показателя Херста при использовании вейвлет-функций с большим числом нулевых моментов. На рис. 4 приведены периодическая составляющая P(t), суммарный сигнал X(t) и спектр вейвлет-энергии суммарного сигнала. В данном случае теоретический показатель H=0,8, величина Ratio = 0.7 . Применена вейв-лет-функция db20 и получена оценка показателя H = 0.808 .
Однако с увеличением значения Ratio и числа периодов во временной реализации требуются вейвлет-функции с все большим числом нулевых моментов, что сложно реализуемо на практике. Кроме того, использование вейвлет-функций, обладающих большим числом нулевых моментов, приводит к менее точному оцениванию показателя Н [16].
1000
0.5
0
-0.5
В работах [11,17] предложен метод корректного оценивания параметра н, если соотношение тренда и фрактального шума является значительным, т.е. значение Ratio=1. В этом случае необходимо предварительно провести анализ спектра вейвлет-энергии сигнала. Если сигнал содержит значительную трендовую или циклическую составляющую, спектр вейвлет-коэффициентов демонстрирует большие значения Ek на соответствующих уровнях разложения k , и при оценивании показателя Н преобладающие в сигнале трендовые составляющие подавляют все остальные.Для этого случая предложено при оценке наклона
графика функции log2(Ej) не учитывать значения трендовой составляющей Ek, тогда оценка показателя Херста будет корректной. На рис. 5 представлены модельный сигнал (Ratio = 4), спектр вейвлет-энергии, полученный с помощью вейвлет-функции db4, и зависимость log2(Ej) с аппроксимирующими прямыми, построенными с учетом и без учета значений Ek при значении k = 8.
Periodic component
0.2 0 -0.2
I
0.5 0 -0.5
I
0.5
0
0
200
400 600
Time series
800
1000
200 400 600 800
Wavelet-energy
1000
0
2
3
4
5
6
7
8
Рис. 4. Периодическая компонента, суммарный сигнал и спектр вейвлет-энергии суммарного сигнала
0 100 200 300 400 500 600 700 800 900 1000 Wavelet-energy
100г
50-
0L
1 2 3 4 5 6 7 8 5I log2 E ' ' ' ' ' '
1
0-5:
2
1 2 3 4 5 6 7 8
E
Рис. 5. Сигнал, спектр вейвлет-энергии, зависимость 1оИ2 (Ер и аппроксимирующие прямые с
учетом (1) и без учета (2) значений Е8
Однако такой подход не приемлем для временных рядов небольшой длины или в случае, когда циклические компоненты попадают сразу в несколько уровней вейвлет-разложения и удаление нескольких значений вейвлет-спектра Ек может привести к неправильному оцениванию показателя Херста. В данной работе предложен метод вейвлет-оценивания пока-
зателя Херста для временных рядов со значительными трендовыми или циклическими составляющими, основанный на использовании пакетного вейвлет-преобразования. Рассмотрим его на примере сигнала, показанного на рис. 6. В данном случае теоретический показатель Херста H=0,8, величина Ratio = 3 . Спектр вейвлет-энергии получен с использованием вейвлет-функции db10. Непосредственное вычисление значения Н приводит к
некорректной завышенной оценке показателя H = 1.04 .
2 1 0 -1 -2
20 15 10 5
0 100 200 300 400 500 600 700 800 900 1000 Wavelet-energy
Time series
0
2
3
4
5
6
7
8
Рис. 6. Сигнал и спектр вейвлет-энергии Очевидно, что основные циклические составляющие попадают в 6-й и 7-й уровни вейвлет-разложения. Детализируем эти уровни, применяя пакетное вейвлет-преобразова-ние с той же вейвлет-функцией db10. Часть дерева пакетного вейвлет-разложения, соответствующая узлам 6-го и 7-го уровней, представлена на рис.7.
Рис.7. Дерево пакетного вейвлет-разложения сигнала, показанного на рис.6 Узел 6.1 дерева пакетного ДВП содержит набор вейвлет-коэффициентов, соответствующих частотному диапазону 6-го уровня разложения обычного ДВП. Узлы 8,4-8,7 пред-
ставляют частотные поддиапазоны с интервалом частот в четыре раза уже диапазона узла 6,1. По формуле (3) рассчитаем значения вейвлет-энергии в узлах 8,4-8,7. Среднее значение этого спектра соответствует значению вейвлет-энергии в узле 6,1, т.е. 6-му уровню разложения ДВП:
E6 = 2 (E7.2 + E7.3 ) = - (E8.3 + E8.4 + E8.5 + E8.6 )•
Аналогичные расчеты выполнены для узла 7,1, т.е. 7-го уровня разложения ДВП. Значения спектра вейвлет-энергии представлены на рис.8. Очевидно, что основная энергия 6-го уровня сосредоточена в узлах 8,5-8,6, а энергия 7-го уровня - в узлах 9,3-9,4. Для наглядности восстановленные компоненты сигнала, которые соответствуют узлам 8,5-8,6, представлены в нижней части рис.8.
Теперь для оценивания показателя Херста в соответствии с отношением (4) в качестве пересчитанных значений энергий 6-го и 7-го уровней возьмем значения
E6 =-(E8.3 + E8.4), E7 =-(E9.5 + E9.6) .
20 10 0
50
200
400
600
800
1000
Рис.8. Часть спектра вейвлет-энергии ДВП, вейвлет-энергия узлов 6,1 и 7,1 пакетного ДВП, восстановленные сигналы узлов 8,5-8,6
На рис. 9 представлен спектр вейвлет-энергии исследуемого сигнала (см. рис. 6)с пересчитанными значениями 6-го и 7-го уровней и полученная зависимость log2(Ej) от номера уровня разложения j. Тангенс угла наклона аппроксимирующей прямой позволяет оценить параметр Н как Н = 0.784, что согласуется со статистическими характеристиками вейвлет-оценок показателя Херста по временным рядам аналогичной длины [16,17].
Wavelet-energy
0.4 0.3 0.2 0.1
1 2 3 4 5 6 7
log2 E
1 2 3 4 5 6 7
6
7
0
0
0
0
Рис. 9. Пересчитанный спектр вейвлет-энергии и соответствующая зависимость l°g2(E j)
Предложенный метод можно применять к любому узлу дерева пакетного ДВП с разной степенью детализации. Таким образом можно получать вейвлет-оценки параметра Херста для самоподобных временных рядов со значительными трендовыми и циклическими составляющими, не теряя значений спектра вейвлет-энергии. Однако надо отметить, что не существует универсального подхода к определению уровня детализации пакетного вейв-лет-дерева. Выбор подходящей вейвлет-функции и количества уровней разложения зависит от временных и частотных характеристик каждого анализируемого сигнала, а также от параметров циклических компонент.
Выводы
Предложен новый метод оценивания показателя Херста, основанный на использовании пакетного вейвлет-преобразования. Применение разработанного метода позволяет определять степень самоподобия для рядов небольшой длины, обладающих существенными трендовыми и циклическими компонентами.
Практическая значимость метода состоит в том, что он основан на алгоритме быстрого вейвлет-преобразования, а это позволяет с минимальными затратами времени и объема информации определять степень самоподобия сложных сигналов, что, в частности, важно для работы в режиме реального времени.
Литература: 1. Федер Е. Фракталы. М.: Мир, 1991. 254 с. 2. Мандельброт Б. Фрактальная геометрия природы. М.: Институт компьютерных исследований, 2002. 656 с. 3. Шелухин О. И. Самоподобие и фракталы. Телекоммуникационные приложения / О. И. Шелухин, А. В. Осин, С. М. Смольский. М. : Физматлит, 2008. 368 с. 4. Шредер М. Фракталы, хаос и степенные ряды. Миниатюры из бесконечного рая / М. Шредер. Ижевск : НИЦ «Регулярная и хаотическая динамика», 2005. 528 с. 5. Clegg R. G. A practical guide to measuring the Hurst parameter / R. G. Clegg // Computing science technical report. 2005. -№ CS-TR-916. 6. Kantelhardt J.W. Fractal and Multifractal Time Series. 2008 [Электронный ресурс]: http:/ /arxiv.org/abs/0804.0747. 7. Малла С. Вэйвлеты в обработке сигналов. М.: Мир. 2005. 671 с. 8. Abry P. Scaling, Fractals and Wavelets / P. Abry, P. Goncalves, J. L. Vehel. - London : John Wiley & Sons, Inc., 2009. 464 p. 9. Abry P. Self-similarity and long-range dependence through the wavelet lens / Р. Abry, P. Flandrin, M. S. Taqqu // Theory and applications of long-range dependence. Birkh^ser, 2003. P. 34-38. 10. Abry P. Wavelet analysis of long-range dependent traffic / P. Abry, D. Veitch // IEEE/ACM Transactions Information Theory. 1998. № 1 (44). Р. 2-15. 11. Дейнеко Ж. В. Оценивание параметра Хёрста для временных рядов с трендом методом вейвлет-преобразования / Ж. В. Дейнеко, Л. О. Кириченко // Системи управлшня, нав^ацп та зв'язку. 2010. Вип. 4 (16). С. 85-89. 12. СмоленцевН.К. Основы теории вейвлетов. Вейвлеты в MATLAB // Н.К. Смоленцев. М.: ДМК Пресс, 2005. 304 с. 13. JeongyH.-D. J. A Comparative Study of Generators of Synthetic Self-Similar Teletrafic / H.-D. J. Jeongy, D. McNickle , K. Pawlikowski // Department of Computer Science and Management, University of Canterbury, 1998. 14. Abry P. The wavelet-based synthesis for the fractional Brownian motion proposed by F. Sellan and Y. Meyer: Remarks and fast implementation / Abry Р., Sellan F. // Appl. and Comp. Harmonic Anal. V. 3(4). 1996. Р. 377-383. 15. MisitiM. Wavelet Toolbox for use with MATLAB / M. Misiti, Y. Misiti, G. Oppenheim, J.-M. Poggi // The MathWorks Inc. 941 p. 2000. 16. Кириченко Л. О. Оценивание самоподобия стохастического временного ряда методом вейвлет-анализа / Л. О. Кириченко, Ж. В. Дейнеко // Радюелектронш i комп'ютерт системи. 2009. № 4 (38). С. 99-105. 17. Kirichenko L. Comparative Analysis for Estimating of the Hurst Exponent for Stationary and Nonstationary Time Series / L. Kirichenko, T. Radivilova, Zh. Deineko // Information Technologies & Knowledge. 2011. Vol. 5, № 4. P. 371-388.
Поступила в редколлегию 12.03.2012 Кириченко Людмила Олеговна, канд. техн. наук, доцент каф. ПМ ХНУРЭ. Научные интересы: вейвлеты, самоподобные и мультифрактальные случайные процессы, хаотическая динамика. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, [email protected]
