П1. Методы сведения к стационарности

Главным этапом в процедуре сведения к стационарности является выделение тренда; его удаление сводится к вычитанию тренда (см. формулы (2.1) и (2.3)) или преобразованию исходной реализации (соотношения (2.2) и (2.4)) к формам (2.1) и (2.3) с последующим вычитанием тренда.

Для того, чтобы успешно произвести операцию выделения и удаления тренда, используют, как правило, три метода. Рассмотрим их.

1. Методы регрессионного анализа^П1 позволяют построить математическую модель, которая наиболее полно описывает функциональный вид тренда. Предполагается, что все реализации изучаемого СП могут быть записаны в форме (2.3):

$$x(t_i)=a(t_i, {\bf b})+\xi (t_i),~~~i=1,\dots , N.$$ (П5)

В этой модели случайного процесса роль зависимой переменной (отклика) играют значения временного ряда $x(t_i)$, а в качестве независимой переменной, влияющей на отклик, выступает дискретное время $t_i$. В формуле (2.5) $x(t_i)$ -- дискретный временной ряд длины $N$, $\xi(t_i)$ -- независимые и одинаково распределенные случайные величины, плотность вероятности которых обычно является гауссовой, $a$ -- функция тренда, ${\bf b}$ -- вектор параметров модели тренда.

Приведем часто применяемые в задачах анализа временных рядов модели трендов $a(t, {\bf b})$.

Наиболее полезной, несмотря на свою простоту, является линейная модель:

$$a(t)=b_0+b_1 t,$$ (П6)

которая описывает, например, изменение характеристик радиофизических устройств после момента включения в течение времени перехода прибора в рабочий режим ("нагревание" устройства).

Для описания нелинейного тренда используется одна из следующих моделей.

Полиномиальная: $a(t)=b_0+b_1t+b_2t^2+...+b_nt^n$; здесь $n$ определяет порядок модели и в практических задачах редко превышает 5. Эта модель описывает плавные неповторяющиеся изменения характеристик процесса.

Логарифмическая: $a(t)=\log (b_0+b_1 t)$. Эта модель применяется, в основном, при анализе переходных процессов.

Логистическая:   $a(t)=b_0/(1+b_1 \exp(-b_2t))$.

Модель Гомперца: $\log(a(t))=b_0-b_1 b_2^t$, $0<b_2<1$.

Логистическая модель и модель Гомперца задают кривые тренда $S$-образной формы и соответствуют процессам, характеристики которых резко возрастают, а затем убывают.

Полигармоническая: $a(t)=\sum_{j,k} b_{jk}\cos(j\omega_kt+\Theta_{jk})$. Эта модель описывает периодические и квазипериодические изменения во временном ряде и применяется в том случае, если известны периоды периодических составляющих процесса, например, из теоретического анализа объекта, генерирующего временной ряд. Если в этой модели параметры $b$ и $\Theta$ рассматривать как случайные равномерно распределенные величины, то она будет описывать поведение стационарного СП. В этом случае можно проводить анализ исходного процесса без удаления полигармонического тренда, что и осуществляется на практике, когда длина временного ряда намного превышает период тренда. Однако, если временной ряд имеет длительность всего несколько периодов тренда, и (или) тренд отличается от чисто периодического процесса, то используется процедура выделения и удаления тренда.

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

$$\sum_{i=1}^{N}(x(t_i)-a(t_i,{\bf b}))^2 \rightarrow min_{\bf b}.$$ (П7)

Рассмотрим процедуру определения параметров тренда на примере линейной модели (2.6). Методом наименьших квадратов найдем оценки параметров $b_0$ и $b_1$, при которых достигается минимум выражения:

$$\sum_{i=1}^{N}(x(t_i)-b_0 -b_1 t_i) ^2 .$$ (П8)

Для определения экстремальных значений (в данном случае минимумов) функции (2.8) приравняем нулю частные производные по $b_0$ и $b_1$ соотношения (2.8). В результате получим алгебраическую систему уравнений относительно неизвестных $b_0$ и $b_1$:

$$\left\{ \begin{array}{c} \sum_{i=1}^{N}[x(t_i)-b_0 - b_1 t_i]=0, \\ \sum_{i=1}^{N}t_i[x(t_i)-b_0 - b_1 t_i]=0. \end{array}\right.$$ (П9)

Решение $(b_0, b_1)$ этой системы с учетом равенства $t_i=iT$ имеет вид:

$$b_0 = \frac{2(2Т-1)\sum_{i=1}^{N}x(i)-6\sum_{i=1}^{N}(i-1)x(i)} {N(N+1)},$$

$$b_1 = \frac{12\sum_{i=1}^{N}(i-1)x(i)-(N-1)\sum_{i=1}^{N}x(i)/2} {TN(N^2-1)}.$$ (П10)

Вычисление коэффициентов $b_0$ и $b_1$ значительно упрощается, если длина ряда $N$ является нечетным числом. В этом случае многие члены выражения (2.10) обращаются в нуль и оно приобретает форму:

$$b_0=\frac{\sum_{i=1}^{N}x(i)}{N},~~~ b_1=\frac{\sum_{i=1}^{N}(i-N/2-1/2)x(i)} {\sum_{i=1}^{N}(1-N/2-1/2)^2}.$$ (П11)

Очевидно, что при четном $N$ последнюю точку временного ряда можно не рассматривать.

Подчеркнем еще раз, что рассмотренный подход применяется к случайным процессам, реализации которых могут быть представлены в виде (2.3). Ясно, что операция удаления тренда сводится в данном случае к вычитанию ряда $a(t_i)$ из исходной последовательности $x(t_i)$, что эквивалентно удалению нестационарного среднего значения.

Допустим теперь, что мы наблюдаем нестационарный процесс вида (2.4): $x(t_i)=a(t_i)u(t_i)$, в котором $u(t_i)$ -- стационарная часть с дисперсией $\sigma ^2_{\xi}$, а $a(t_i)$ -- некоторая положительная неслучайная последовательность. В этом случае дисперсия процесса изменяется во времени: $\sigma_x^2 (t)= \sigma ^2_{\xi}a^2(t)$. Нестационарность, как видно, проявляется в нерегулярном изменении энергии (дисперсии) ряда, поэтому перед применением описанных выше моделей необходимо преобразовать временной ряд к форме (2.3), т.е. изменить шкалу значений ряда, в которой они измерены. Чаще других используется переход к логарифмической шкале:

$$y(t_i)=\log(x(t_i)+c);$$ (П12)

здесь $c$ -- некоторая  константа,   которая  выбирается  из  условия: $(x(t)+c)>0$ для всех $t_i$. Переход к логарифмической шкале приводит к формуле:

$$\begin{eqnarray*} y(t_i)=\log(a(t_i))+\log(u(t_i)), \end{eqnarray*}$$

в которой $\log(a(t_i))$ и $\log(u(t_i))$ являются нестационарной и стационарной частями, соответственно. Таким образом, мы получили модель тренда вида (2.3), параметры которой могут быть найдены методом наименьших квадратов.

2. Метод скользящих средних. Если априорная информация о характере тренда отсутствует, то для его удаления используют метод скользящих средних. Этот метод основан на представлении нестационарной части временного ряда $a(t_i)$ в виде последовательности средних значений исходного ряда, вычисленных на коротком временном интервале, центр которого "скользит" вдоль всего ряда. По сути дела, проводится процедура усреднения последовательности $x(t_i)$ в плавающем окне, в результате чего ряд скользящих средних $a(t_i)=\bar{x}(t_i)$, являющийся трендом, ведет себя более гладко, чем исходный. Далее вычитанием $\bar{x}(t_i)$ из исходного ряда осуществляется переход к стационарной последовательности: $u(t_i)=x(t_i)-\bar{x}(t_i)$.

Рассмотрим процедуру вычисления cкользящих средних дискретной реализации $x(t_i)$. Выберем длину окна усреднения $p=2m+1$. Скользящие средние в этом случае можно вычислить по формуле:

$$a(t_i)=\bar{x}(t_i)=\frac{1}{2m+1} (x(t_{i-m})+...+x(t_{i-1})+x(t_i)+x(t_{i+1})+...+x(t_{i+m})).$$ (П13)

Ряд $\bar{x}(t_i)$ содержит на $p$ меньше членов, чем исходный, что приводит к появлению неопределенности в поведении скользящих средних на краях реализации. Поэтому операция удаления тренда выполняется для более короткого (на $p$ членов) временного ряда $x(t_i)$; стационарная последовательность $u(t_i)$ при этом имеет длину $N-p$.

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

3. Метод разностных операторов. Этот метод заключается в переходе от исходной временной последовательности к ряду разностей соседних значений:

$$y(t_i)=x(t_i)-x(t_{i-1})=\nabla x(t_i).$$ (П14)

Выражение (2.14) называется разностным оператором первого порядка.

Разностный оператор второго порядка имеет вид:

$$y(t_i)=\nabla^2 x(t_i)=x(t_i)-2x(t_{i-1})+x(t_{i-2}).$$ (П15)

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

Рассмотрим действие метода разностных операторов на временной ряд, имеющий линейный тренд и интервал дискретизации, равный единице: $t_i-t_{i-1}=1$:

$$y(t_i) = \nabla x(t_i)=x(t_i)-x(t_{i-1})=$$ (П16)

$$~ = b_0+b_1t_i+u(t_i)-b_0-b_1t_{i-1}-u(t_{i-1})=b_1+u(t_i)-u(t_{i-1}).$$

В отличие от исходной последовательности $x(t_i)$, преобразованный ряд $y(t_i)$ короче на один член и не содержит нестационарности. В результате такого преобразования структура стационарной компоненты изменяется.

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