Любушин А.А. Фрактальный анализ временных рядов

Подождите немного. Документ загружается.

f) Значение параметра Херста 1 для самоподобного процесса. Действительно, имеем:

H ≤

|(2)| |(2) (1) (1)| |(2) (1)| |(1)|2 |(1)|

MZ MZ Z Z MZ Z MZ MZ=−+≤−+= (37)

С другой стороны, в силу самоподобия |(2)|2 |(1)| 2 2 1

MZ MZ H

⇒≤⇒≤ .

g) Если , то

1H =

() () (,)

ZsZt st st

=Γ = . Отсюда

22 22 22 2

(() (1)) ()2 (()(1)) (1) ( 2 ) 0MZt tZ MZ t tMZtZ tMZ t tt t

−= − + =−+= (38)

то есть () (1)

ttZ= почти всюду.

Пусть ()

t Hsssi∈ , и ()01H<<

t - гауссовский процесс. Тогда ()

t называется

фрактальным броуновским движением и обозначается ()

t . Если , то 0.5H =

()

является обычным броуновским движением или винеровским процессом. Нетрудно

получить, что в случае, если

0.5

(,) min(| |,| |)st s t

Γ=

() ()sign t sign s

и , если

. Одно из возможных представление fBm:

0.5

(,) 0stΓ=

() ()sign t sign s≠

0.5 0.5

0.5

() [( ) ( ) ] ( ); max(0, )

tk tu u dBuu u

+∞

−−

++ +

−∞

=−−− =

∫

(39)

Здесь

k - нормировочная постоянная. Если выбрать:

2(1.5)

(0.5)(22)

() , 0,(0.5) ,(1)1,(1)

uetdtu n n

∞

−−

⋅Γ −

Γ+ ⋅Γ−

Γ= > Γ = Γ= Γ+=

∫

(40)

то . Хаусдорфова размерность реализаций

() 1

MB t

= ()

t равна 2.

DH=−

Рассмотрим приращения процесса ()

t : ( 1) ( )

kH H

Bk Bk

+− . Эта случайная

последовательность называется

фрактальным гауссовским шумом. Для фрактального шума

справедливы утверждения:

(i)

стационарна.

(ii)

MX MX

(iii)

22 2 2

() (| 1| 2| | | 1| )/2

HH H

XX i i k

kMXX k k k

γσ

==+−+−

(iv) Пусть .

0k ≠

Тогда:

0.5 ( ) 0;

00.5 ()

0.5 1 ( ) 0

=⇒

<< ⇒ <

<<⇒ >

(v) Пусть . Тогда 0.5H ≠

222

() (2 1)| | , | |

kHHk k

γσ

−

⋅→∼ ∞

(vi) Пусть

() ()

XX XX

Ske

ωγ

+∞

−

=−∞

∑

- спектральная плотность стационарной случайной

последовательности

; ( ) ( )

XX XX

kSed

−

∫

. Тогда , т.е. если

(2 1)

() , 0

ωω ω

−−

→∼

0.5 1 ( ) , 0

<<⇒ →∞→ и сигнал

носит низкочастотный характер.

Соответственно, спектральная оценка конечной выборки из фрактального броуновского

движения будет иметь поведение:

(2 1)

() , 0

ωω ω

−+

→∼

Оценка постоянной Херста по экспериментальным данным. Наиболее очевидная

оценка непосредственно из определения показателя, данного самим автором. Пусть

( ), 1,...,

tt N – анализируемый временной ряд;

– длина скользящего временного

окна;

- номер отсчета правого конца скользящего окна, то есть мы рассматриваем

моменты времени , которые удовлетворяют условию

t 1

−

+≤≤Lt

. Пусть - длина

внутреннего временного окна, которое используется внутри текущего основного окна для

операций усреднения. Мы рассматриваем длины внутренних окон, удовлетворяющие

условию: . Пусть

/5≤sL

()

(

+−

∑

xxut

(41)

− выборочная оценка среднего значения на интервале длиной отсчетов, который лежит

внутри текущего основного окна и начинается в точке . Следующий шаг состоит в

вычислении отклонений от среднего значения (41), их накопленной суммы и размаха

накопленной суммы:

() () () () () () ()

,,,,,,

() () , () (), max () min ()

Δ=− =Δ = −

∑

su su su su su su su

txtx t xvR t t

ττττττ

ξξ

(42)

для [, 1], [ 1, 1]

∈+− ∈−+−+tuus u L s

. Далее оцениваются дисперсия и среднее значение

отношения размаха к стандартному отклонению:

()

() 2 () 2 ()

()

() ( ()), ()

(1)

−+

=Δ =

−+

∑∑

−+

su su

xt RSs

sLs

τττ

(43)

Показатель Херста ()

в текущем временном окне оценивается как наклон кривой

прямой линейной регрессии между значениями

()

ln( ( ))

и . Оценка может получена

как в скользящем временном окне, так и по всей выборке. Заметим, что если ставится задача

оценки постоянной Херста для процесса

ln( )s

()

t Hsssi

∈

, то все вышеизложенные операции

следует производить после перехода к его приращениям, то есть для

() ( 1) ()

=+−

tZt Zt.

Другие оценки могут быть основаны на свойстве роста спектра мощности по степенному

закону: при уменьшении частоты. Если использовать не обычные

спектральные оценки, а разложения по системе ортогональных нормированных компактных

базисных функций – вейвлетам [1], то для процесса

(2 1)

() , 0

−+

→∼

ωω ω

()

t Hsssi

∈

оценку постоянной Херста

можно найти по скорости роста средних значений квадратов модулей вейвлет-

коэффициентов:

()

()2 ()

||/

∑

WcN

(44)

Здесь

()

- коэффициенты ортогонального дискретного вейвлет-разложения выборки

самоподобного временного ряда,

1,...,

= m

- номер уровня детальности разложения,

()

число вейвлет-коэффициентов на уровне детальности

() ( )

−

≤

. Тогда, аналогично

соотношению для скорости роста спектра мощности,

()

∼

αα

, где s

- характерный

временной масштаб уровня детальности

. Поскольку

(1)

=÷s

αα

, то отсюда следует, что

log ( )

∼

(45)

Таким образом, значение коэффициента наклона прямой, подогнанной методом

наименьших квадратов к парам значений

(log ( ), )W

, дает оценку для величины 21

H .

Оценка спектра сингулярности по экспериментальным данным. В настоящее время

существуют 2 подхода для оценки спектров сингулярности временного ряда. Первый метод

появился раньше и основан на анализе цепей точек максимума модулей непрерывных

вейвлет-преобразований с вейвлетами, обычно равными производной той или иной степени

от функции плотности распределения Гаусса [1]. Второй подход более близок к технике

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

ее длины. В последнее время был разработан и активно применяется в различных

приложениях метод анализа флуктуаций после исключения масштабно-зависимых трендов –

Detrended Fluctuation Analysis (DFA) [4]. Сравнительный опыт применения методов

показывает, что метод DFA является более надежным и устойчивым. В то же время для

специального вида самоподобных сигналов, которые могут содержать плато постоянных

значений (типа известной «чертовой лестницы», конструируемой на базе канторовского

множества), метод DFA неприменим и оценка, основанная на непрерывных вейвлет-

преобразованиях, имеет преимущества. Ниже будет использован только DFA и приведены

основные конструкции метода.

Пусть ()

t – случайный процесс. Определим в качестве меры (, )

изменчивости

сигнала ()

t на интервале [, ]tt

+ модуль его приращения: ( , ) | ( ) ( )|=+−

txt xt

δδ

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

q :

(,) {( (,))}

MqM t

δμδ

(46)

Случайный процесс называется масштабно-инвариантным, если при

()

(,) | |

δδ

∼

→

, то есть существует предел:

ln ( , )

() lim

ln | |

→

= (47)

Заметим, что в определении (46-47) величина меры (, )

может быть взята также как

размах, что ближе к традиционным конструкциям Херста:

(, ) max ( ) min ( )

≤≤+

−

tut

txu

(48)

Если зависимость

()q

является линейной:

()qHq

, где , то

процесс называется монофрактальным. В частности, для классического броуновского

движения . Возведение в степень подчеркивает различные типы поведения

,0 1H const H=<<

0.5H =

сигнала: если 0, то в значение меры (46) основной вклад вносят интервалы времени с

большими отклонениями от тренда, а если 0

q , то интервалы времени с малыми

вариациями.

Для вычисления функции ()

по конечной выборке из временного ряда ( ), 1,...,

tt N

можно применить метод DFA. Пусть - число отсчетов, ассоциированное с варьируемым

масштабом

=Δ. Разобьем выборку на непересекающиеся малые интервалы длиной

отсчетов:

()

{ :1 ( 1) , 1,...,[ / ]}

tkstksk Ns=+−≤≤ = (49)

и пусть

()

( ) (( 1) ), 1,...,=−+ =

txk stt s (50)

участок временного ряда ()

t , соответствующий интервалу

()

. Пусть - полином

порядка , подогнанный методом наименьших квадратов к сигналу

(, )

()

t . Рассмотрим

отклонения от локального тренда:

(, ) () (, )

( ) ( ) ( ), 1,...,

sm s sm

kkk

ytytpttΔ=− =s (51)

и вычислим значение:

[/]

() (,) (,)

(,) (max () min ()) [ /]

msmsmq

Zqs y t y t Ns

≤≤

⎛⎞

=Δ−Δ

⎜⎟

⎝⎠

∑

(52)

которое будем рассматривать как оценку для . Процедура устранения тренда на

каждом малом участке длиной отсчетов необходима в случае наличия в сигнале трендов

внешнего происхождения. Определим теперь функцию как коэффициент линейной

регрессии между значениями

((,))

()hq

()

ln( ( , ))

qs и : ln( )s

() ()

(,)

mhq

qs s∼ . Очевидно, что

() ()qqhq

= , а для монофрактального процесса ()h q H const

= .

Следующим шагом в мультифрактальном анализе после определения функции

()q

является вычисление спектра сингулярности

()F

, который является с фрактальной

размерности множества точек, в окрестности которых показатель Гельдера-Липшица для

случайных реализаций процесса

()

t равен

, то есть таких точек , для которых t

|( ) ()|| |, 0

+− →∼xt xt

δδδ

. Стандартный подход состоит в вычислении статистической

суммы Гиббса:

[/]

(, ) (, )

(,) (max () min ())

msm

Wqs y t y t

≤≤

=Δ−Δ

∑

(53)

и определения показателя массы ()

из условия , после чего спектр

()

(,)

Wqs s

∼ ()F

вычисляется согласно формуле:

() max{min( ()),0}

Fqq

ατ

=− (54)

Сравнивая (52) и (53), нетрудно заметить, что () () 1 () 1

qqqhq

−= −. Таким образом,

( ) max { min( ( ( )) 1, 0}

Fqhq

=−+.

Для монофрактального процесса, когда

()h q H const

= , получаем, что ()1FH

() 0FH

=∀≠. В частности, положение и ширина носителя спектра ()F

, то есть

значения

min max max min

αααα

Δ= − и

- то значение, которое доставляет функции ()F

максимум:

()max()F

, являются характеристиками шума. Величину

можно

назвать обобщенным показателем Херста. Для монофрактального сигнала теоретически

значение

Δ должно быть равно нулю, а

, но на практике, вследствие конечности

выборки, такие условия не выполняются. Что же касается значения

()F

, то оно равно

фрактальной размерности точек, для окрестности которых выполняется масштабирующее

соотношение .

()

(,) | |

δδ

∼

Если оценивать спектр ()

в скользящем временном окне, то его эволюция может дать

информацию об изменении структуры хаотических пульсаций ряда. Обычно

()1F

, но

встречаются окна, для которых

()1F

. Напомним, что в общем случае (не только для

анализа временных рядов) величина

()F

равна фрактальной размерности носителя

мультифрактальной меры.

Пример анализа данных. На рис.1(а) представлен график временного ряда наблюдений

за вариациями электротеллурических потенциалов на Камчатке из работы. Общая

продолжительность наблюдений составляет 4 года 8 месяцев (с 01.10.1996 по 23.06.2001),

интервал взятия отсчетов 1 час, общее число отсчетов N=41447. Рис.1(б) представляет

приращения временного ряда на рис.1(а).

Рис.1. Графики исходного ряда наблюдений и его приращений.

Рис.2 иллюстрирует три различных метода оценки постоянной Херста для временного

ряда, представленного на рис.1. На рис.2(а) оценка получена как коэффициент наклона

прямой наилучшего приближения (штриховая линия) к графику логарифма оценки спектра

мощности в зависимости от логарифма периода. Оценка спектра мощности получена

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

скользящем временном окне длиной 2048 отсчетов с использованием модели авторегрессии

200-го порядка внутри каждого окна. На графике рис.2(а) видны спектральные пики,

соответствующие 8, 12 и 24-часовым периодическим компонентам. На рис.2(б) использована

оценка (44)-(45), ортогональный вейвлет выбран из условия минимума энтропии

распределения квадратов вейвлет-коэффициентов [1] и им оказался вейвлет Добеши 18-го

порядка (обнуляющий 9 первых моментов). Наконец, рис.2(б) иллюстрирует применение

статистики (43).

Рис.2. Оценка показателя Херста тремя различными методами.

Эти методы дали значения постоянной Херста 0.37, 0.40 и 0.46. Различие в оценках не

очень существенно, но не пренебрежимо мало. Оценка на рис.2(а) представляется наименее

надежной из-за влияния монохроматических пиков. Кроме того, следует учесть, что сигнал,

возможно, является мультифракталом – в этом случае отличия возникают просто вследствие

различия методов и различий в «реакции» на них данных, имеющих на самом деле иную

природу. Для выяснения этого вопроса оценим спектр сингулярности для приращений

временного ряда (рис.1(б)) – см. рис.3.

Рис.3. Оценка показателя спектра сингулярности (а) и зависимостей

() ()

(,)

mhq

qs s∼ для

различных показателей степени

q (б)

На рис.3(б) представлены графики как самих величин

()

lg( ( , ))

qs , так и линейных

трендов, подогнанных к ним для 9 различных значений степени . Наклон графиков

линейных трендов есть ничто иное, как функция (), по значениям которой вычисляется

показатель массы

() () 1

=qqhq−

, а далее, согласно формуле (54), и сам спектр

сингулярности. Рис.3(б) демонстрирует хорошее соответствие линейному закону

зависимости

()

lg( ( , ))

qs от , то есть, самоподобность сигнала. На рис.3(б) приведены

графики зависимостей с подогнанными линейными трендами для последовательных

значений степени (сверху вниз): 10, 7.5, 5, 2.5, -0.05, -2.5, -5, -7.5 и -10. Масштаб

изменялся начиная с минимального значения 30 отсчетов (30 часов) – чтобы усреднить

периодические компоненты, максимальный период которых равен 24 часам. Следовательно,

как это видно из рис.3(б), сигнал действительно является мультифракталом.

lg( )s

Для получения временных рядов эволюции значений

было выбрано скользящее окно

длиной 672 часовых отсчетов, то есть 28 суток, взаимное смещение соседних окон равнялось

24 отсчетам или 1 суткам. Масштабно-зависимые тренды не устранялись, так как

рассматривался ряд в приращениях. Функция () в зависимости

() ()

(,)

mhq

qs s∼

оценивалась в каждом окне для масштабов , изменяющихся от минимального значения 30

отсчетов до максимального, равного одной пятой длины окна. При длине 672 отсчетов,

таким образом, максимальный масштаб равен 137 часов. Кроме того, для этих же временных

окон были получены оценки на основе использования статистик (43) для ряда в приращениях

и (44)-(45) и для исходного ряда. При этом, при построении оценки на основе ортогональных

вейвлет-разложений производилось предварительное устранение тренда внутри каждого

окна полиномом 4-го порядка. Эта предварительная операция удаляла статистически

незначимые низкочастотные вариации и, что более важно, уменьшало искажения значений

вейвлет-коэффициентов на концах выборки вследствие ее конечности. В каждом окне

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

квадратов вейвлет-коэффициентов.

Из рис.4(а) и (б) видно, что оба метода, (43) и (44)-(45) дают, фактически одно и то же, за

исключением взаимного смещения, которое может быть вызвано несоответствием данных

монофрактальному поведению.

Что же касается прочих графиков на рис.4, то они могут быть использованы при

мониторинге объектов для выделения аномальных признаков. Первым признаком

аномального поведения может служить «провал» вниз значения

()F

- это может

свидетельствовать о появлении во временном ряду составляющей, поведение которой

достаточно сильно отличается от поведения случайной самоподобной кривой. Кроме того,

вариации наиболее типичного показателя Гельдера-Липшица

и ширины носителя

спектра сингулярности могут также нести информацию об изменении свойств ряда, которые

скрыты от их выделения более простыми средствами анализа типа линейной фильтрации.

Рис.4. (а) – оценки вариаций постоянных Херста в скользящем временном окне длиной 672

часа RS-методом (формула (43)); (б) – оценка с использованием ортогональных вейвлет

разложений содержимого каждого из окон (формулы (44) и (45)). Прочие графики

представляют изменения параметров спектра сингулярности (обозначения – в тексте).