Санкт-Петербургский государственный университет
Математико-механический факультет
Вычислительная стохастика и статистические модели
Круглова Валентина Викторовна
Идентификация кратных корней знаменателя спектральной
плотности стационарных моделей
Бакалаврская работа
Научный руководитель:
к. ф.-м. н., доцент Т. М. Товстик
Рецензент:
к. ф.-м. н., доцент Н. М. Москалева
Санкт-Петербург
2016
Saint Petersburg State University
Applied Mathematics and Computer Science
Computational Stochastics and Statistical Models
Kruglova Valentina Victorovna
Identification of multiple roots for denominator in spectral
density of stationary models
Bachelor’s Thesis
Scientific Supervisor:
Associate Professor T. M. Tovstik
Reviewer:
Associate Professor N. M. Moskaleva
Saint Petersburg
2016
3
Содержание
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Глава 1.
4
Идентификация кратных корней знаменателя
спектральной плотности стационарных процессов . . . . . . . . . . . . .
6
1.1.
Основные определения . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.2.
Оценка коэффициентов стохастического уравнения . . . . . . . . . . . .
10
1.3.
Определение порядка процесса авторегрессии . . . . . . . . . . . . . . . .
13
1.4.
Выявление кратных корней характеристического полинома . . . . . . . .
14
1.5.
Оценка погрешностей параметров 𝛼
^и𝜔
^ у процесса АР(2) . . . . . . . .
17
1.6.
Моделирование процесса с дискретным параметром . . . . . . . . . . . .
20
Глава 2.
Практическое применение . . . . . . . . . . . . . . . . . . . . . . .
22
2.1.
Стационарный процесс с дискретным параметром . . . . . . . . . . . . .
22
2.2.
Подбор стационарной модели по методу Юла–Уокера . . . . . . . . . . .
26
2.3.
Подбор стационарной модели с помощью пошагового метода . . . . . . .
34
2.4.
Сравнительная характеристика методов . . . . . . . . . . . . . . . . . . .
42
2.5.
Приложение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
4
Введение
В современном мире статистическая обработка временных рядов используется для
выявления изменения тенденций, периодических закономерностей, для прогноза ряда
или его моделирования. В данной работе основное внимание было уделено исследованию
стационарных случайных процессов с дробно–рациональной спектральной плотностью,
а именно, изучению методов выявления кратных корней знаменателя спектральной
плотности стационарных моделей.
Определение. Случайный процесс 𝑋(𝑡) называется стационарным в широком смыс
ле, если математическое ожидание процесса является константой, то есть E𝑋(𝑡) =
𝑐, кроме того, E|𝑋(𝑡)|2 < ∞ и корреляционная функция 𝑅(𝑡, 𝑠) = E𝑋(𝑡)𝑋(𝑠) зависит
лишь от разности 𝑡 − 𝑠: 𝑅(𝑡, 𝑠) = 𝑅(𝑡 − 𝑠).
На практике очень часто встречаются случайные процессы, квадратичные харак
теристики которых остаются неизменными во времени. Именно эти процессы и назы
ваются стационарными. Исследование подобных процессов осуществляется с помощью
методов теории вероятностей и математической статистики.
Основная цель работы состоит в том, чтобы исследовать возможность идентифи
кации и оценки кратности корней знаменателя спектральной плотности стационарных
моделей. Поставленный вопрос имеет важную роль, поскольку наличие кратного корня
говорит о том, что в данной точке существует пик у спектральной плотности, то есть
происходит резонанс.
Изучение свойств спектральной плотности стационарных процессов дает возмож
ность характеризовать спектральный (или частотный) состав процесса. Этот факт ши
роко используется при изучении различных сигналов и процессов.
В ходе работы были поставлены следующие задачи:
∙ Моделирование процесса с дискретным параметром
∙ Статистическая обработка полученных результатов
∙ Выявление возможности оценки кратности корней знаменателя спектральной плот
ности построенных моделей.
5
В ходе работы были исследованы методы оценивания параметров различных про
цессов, а также возможность идентификации кратных корней знаменателя спектраль
ной плотности построенных моделей. Теоретический материал по данному вопросу пред
ставлен в главе 1, а результаты статистической обработки описаны в разделах 2.2 и 2.3.
Кроме того, были исследованы общие принципы моделирования, применимые к
процессам из рассматриваемого класса. Способы моделирования случайных процессов с
дробно–рациональными спектральными плотностями подробно изложены в книге [1]. В
разделе 1.6 приведены некоторые изученные методы моделирования, которые подробно
представлены в книге [2]. Данный теоретический материал был реализован на практике
с помощью языка программирования R, широко используемого для статистического
анализа, и полученные результаты приведены в разделе 2.1.
6
Глава 1
Идентификация кратных корней знаменателя
спектральной плотности стационарных процессов
1.1. Основные определения
Сформулируем некоторые определения и теоремы, играющие ведущую роль на
протяжении всей работы.
Определение 1.1. Случайный процесс 𝑋(𝑡), 𝑡 ∈ 𝑇 называется процессом с непрерыв
ным параметром, если 𝑇 является интервалом на R. Если же 𝑇 — конечное или
счетное множество, то процесс 𝑋𝑡 = 𝑋(𝑡), 𝑡 = 0, ±1, ±2, . . . , называется процессом с
дискретным параметром.
Определение 1.2. Процесс 𝑋𝑡 называется процессом авторегрессии с остатками
скользящего среднего (APCC(𝑛, 𝑚)), если он удовлетворяет разностному стохасти
ческому уравнению
𝑋𝑡 + 𝑞1 𝑋𝑡−1 + · · · + 𝑞𝑛 𝑋𝑡−𝑛 = 𝑝0 𝜉𝑡 + 𝑝1 𝜉𝑡−1 + · · · + 𝑝𝑚 𝜉𝑡−𝑚 ,
(1.1)
где 𝜉𝑘 - независимые случайные величины с E𝜉𝑘 = 0, E𝜉𝑘 𝜉𝑗 = 𝛿𝑘𝑗 , а 𝛿𝑘𝑗 — символ Кроне
кера. В частности, если 𝑚 = 0, то есть
𝑋𝑡 + 𝑞1 𝑋𝑡−1 + · · · + 𝑞𝑛 𝑋𝑡−𝑛 = 𝑝0 𝜉𝑡 ,
(1.2)
то процесс 𝑋𝑡 называется процессом авторегрессии 𝑛–го порядка (AP(𝑛)). Если же
𝑛 = 0, то 𝑋𝑡 удовлетворяет уравнению
𝑋𝑡 = 𝑝0 𝜉𝑡 + 𝑝1 𝜉𝑡−1 + · · · + 𝑝𝑚 𝜉𝑡−𝑚
(1.3)
и называется процессом скользящего среднего 𝑚–го порядка (CC(𝑚)).
Теорема 1.1 ([2]). Если все корни характеристического полинома 𝑞(𝑧) =
∑︀𝑛
𝑘=0 𝑞𝑘 𝑧
𝑘
,
𝑞0 = 1, по модулю больше единицы, то процесс 𝑋𝑡 АРСС(𝑛, 𝑚) стационарен и имеет
спектральную плотность следующего вида
⃒
⃒2
1 ⃒⃒ 𝑝(𝑒−𝑖𝜆 ) ⃒⃒
𝑓 (𝜆) =
,
2𝜋 ⃒ 𝑞(𝑒−𝑖𝜆 ) ⃒
(1.4)
7
где полином 𝑝(𝑧) =
∑︀𝑚
𝑘=0
𝑝𝑘 𝑧 𝑘 имеет действительные коэффициенты и корни вне еди
ничного круга.
Спектральная плотность процесса с непрерывным параметром имеет вид
⃒
⃒2
1 ⃒⃒ 𝑃 (𝑖𝜆) ⃒⃒
.
𝑓 (𝜆) =
2𝜋 ⃒ 𝑄(𝑖𝜆) ⃒
(1.5)
При этом полиномы
𝑃 (𝑧) =
𝑚
∑︁
𝑏𝑘 𝑧 𝑚−𝑘 ,
𝑚 < 𝑛,
𝑎𝑘 𝑧 𝑛−𝑘 ,
𝑎0 = 1.
𝑘=0
𝑄(𝑧) =
𝑛
∑︁
𝑘=0
имеют действительные коэффициенты, а их корни лежат в левой полуплоскости.
Для процессов со спектральной плотностью вида 1.4 и 1.5 корреляционная функ
ция имеет следующий вид
𝑅(𝑡) =
𝑠
∑︁
𝑒
−𝛼𝑘 |𝑡|
𝑘=1
𝑣𝑘
∑︁
𝑗
𝑗
(𝐴𝑘𝑗 |𝑡| cos(𝜔𝑘 𝑡) + 𝐵𝑘𝑗 |𝑡| sin(𝜔𝑘 |𝑡|)) +
𝑗=0
𝑝
∑︁
𝑒
−𝛽𝑘 |𝑡|
𝑘=1
𝑤𝑘
∑︁
𝐶𝑘𝑗 |𝑡|𝑗 .
(1.6)
𝑗=0
Заметим, что в формуле (1.6) степень у |𝑡| соответствует кратности корня в знаме
нателе факторизованной спектральной плотности.
Корреляционная функция (1.6) и спектральная плотность процессов с непрерыв
ным и дискретным параметрами связаны соответствующими соотношениями
1
𝑓 (𝜆) =
2𝜋
+∞
Z
𝑒−𝑖𝜆𝑡 𝑅(𝑡)𝑑𝑡, 𝑓 (𝜆) =
−∞
∞
Z
Z𝜋
𝑒𝑖𝜆𝑡 𝑓 (𝜆)𝑑𝜆,
𝑅(𝑡) =
−∞
+∞
1 ∑︁
𝑅𝑗 𝑒−𝑖𝜆𝑗 ;
2𝜋 𝑗=−∞
𝑒𝑖𝜆𝑡 𝑓 (𝜆)𝑑𝜆.
𝑅𝑗 =
−𝜋
Пример 1.1. В данной работе рассматривается стационарный процесс АР(4), имею
щий в знаменателе спектральной плотности два комплексно–сопряженных корня вто
рой кратности. Корреляционная функция и спектральная плотность для этого процесса
имеют вид
𝑅(𝑡) = 𝑒−𝛼|𝑡|
(︁(︀
)︀ (︀
)︀ )︁
𝐴0 cos(𝜔𝑡) + 𝐵0 sin(𝜔|𝑡|) + 𝐴1 cos(𝜔𝑡) + 𝐵1 sin(𝜔|𝑡|) |𝑡|
(︂
)︂
𝑈0 + 𝑈1 cos 𝜆 + 𝑈2 cos(2𝜆) + 𝑈3 cos(3𝜆)
1
𝑓 (𝜆) =
,
2𝜋 𝑉0 + 4𝑉1 cos(𝜆) + 4𝑉2 cos(2𝜆) + 4𝑉3 cos(3𝜆) + 2𝑉4 cos(4𝜆)
(1.7)
(1.8)
8
где
)︀
)︀
(︀
(︀
𝑈0 = 𝐴0 (1 − 𝑒−4𝛼 ) 1 + 4𝑒−2𝛼 + 𝑒−4𝛼 + 4𝐴1 1 − 𝑒−4𝛼 + 2𝑒−6𝛼 +
)︀
(︀
+ 4𝑒−2𝛼 𝐴0 (1 − 𝑒−4𝛼 ) + 2𝐴1 𝑒−4𝛼 cos(2𝜔) +
)︀
(︀
+ 4𝑒−2𝛼 (1 + 𝑒−2𝛼 ) −𝐵0 (1 + 𝑒−2𝛼 ) + 4𝐵1 𝑒−2𝛼 sin(2𝜔) − 2𝐵0 𝑒−4𝛼 sin(4𝜔);
(︂ (︂
)︂
−2𝛼
−2𝛼
−4𝛼
−2𝛼
−2𝛼
−𝛼
𝐴0 (1 − 𝑒 )(3 + 4𝑒
+ 3𝑒 ) cos(𝜔) − 𝑒 (1 − 𝑒 ) cos(3𝜔) +
𝑈1 = 2𝑒
(︂
)︂
−2𝛼
−4𝛼
−6𝛼
−2𝛼
−2𝛼
+ 𝐴1 − (5 − 7𝑒
−𝑒
+ 6𝑒 ) cos(𝜔) + 𝑒 (3 − 5𝑒 ) cos(3𝜔) +
(︂
+ 𝐵0 2𝑒−2𝛼 (1 + 𝑒−2𝛼 ) cos(𝜔) + (1 + 𝑒−2𝛼 )3 sin(𝜔) + 𝑒−2𝛼 (1 + 𝑒−2𝛼 ) sin(3𝜔) +
)︂
(︂
−2𝛼
−2𝛼
+ 2𝑒 (1 + 𝑒 ) cos(3𝜔) + 𝐵1 − (1 + 3𝑒−2𝛼 + 5𝑒−4𝛼 + 7𝑒−6𝛼 ) sin(𝜔) +
)︂)︂
−2𝛼
−2𝛼
+ 𝑒 (3 + 5𝑒 ) sin(3𝜔) ;
(︂
(︂
)︂
−2𝛼
−4𝛼
−4𝛼
−4𝛼
𝑈2 = 2𝑒
𝐴0 (1 − 𝑒 )(2 + cos(2𝜔)) + 𝐴1 4𝑒
− 2(1 − 3𝑒 ) cos(2𝜔) −
)︂
−2𝛼
−4𝛼
−2𝛼
−2𝛼
− 𝐵0 (1 + 4𝑒
+ 𝑒 ) sin(2𝜔) + 2𝐵1 (1 + 𝑒 )(1 + 3𝑒 ) sin(2𝜔) ;
(︂
−3𝛼
𝑈3 = 2𝑒
𝐴0 (1 − 𝑒−2𝛼 ) cos(𝜔) + 𝐴1 (3 − 5𝑒−2𝛼 ) cos(𝜔) + 𝐵0 (1 + 𝑒−2𝛼 ) sin(𝜔) −
)︂
−2𝛼
− 𝐵1 (3 + 5𝑒 ) sin(𝜔) ;
(︀
)︀2
(︀
)︀2
𝑉0 = 1 + 4𝑒−2𝛼 + 𝑒−4𝛼 + 8𝑒−2𝛼 1 + 𝑒−2𝛼 cos(2𝜔) + 2𝑒−4𝛼 cos(4𝜔);
)︂
(︂
(︀
)︀
−2𝛼
−4𝛼
−𝛼
−2𝛼
−2𝛼
+𝑒
cos(𝜔) + 𝑒 cos(3𝜔) ;
𝑉1 = −2𝑒 (1 + 𝑒 ) 1 + 4𝑒
(︂
)︂
(︀
)︀
−2𝛼
−2𝛼 2
−2𝛼
−4𝛼
𝑉2 = 𝑒
2(1 + 𝑒 ) + 1 + 4𝑒
+𝑒
cos(2𝜔) ;
𝑉3 = −2𝑒−3𝛼 (1 + 𝑒−2𝛼 ) cos(𝜔);
𝑉4 = 𝑒−4𝛼 .
Дисперсия рассматриваемого процесса равна коэффициенту 𝐴0 . Рассмотрим слу
чай, когда дисперсия процесса равна 1, то есть 𝐴0 = 1.
Поскольку рассматривается процесс АР(4), то в числителе 𝑓 (𝜆) должна быть кон
станта, не зависящая от 𝜆. По этой причине все три коэффициента 𝑈1 , 𝑈2 , 𝑈3 приравни
ваются к нулю, что приводит к следующей системе линейных алгебраических уравнений
относительно 𝐴1 , 𝐵0 , 𝐵1 :
9
⎧ (︂
)︂
⎪
−2𝛼
−4𝛼
−6𝛼
−2𝛼
−2𝛼
⎪
𝐴1 − (5 − 7𝑒
−𝑒
+ 6𝑒 ) cos(𝜔) + 𝑒 (3 − 5𝑒 ) cos(3𝜔) +
⎪
⎪
⎪
⎪
(︂
⎪
⎪
⎪
⎪
⎪
+𝐵0 2𝑒−2𝛼 (1 + 𝑒−2𝛼 ) cos(𝜔) + (1 + 𝑒−2𝛼 )3 sin(𝜔) + 𝑒−2𝛼 (1 + 𝑒−2𝛼 ) sin(3𝜔)+
⎪
⎪
⎪
)︂
(︂
⎪
⎪
⎪
⎪
−2𝛼
−2𝛼
⎪
+2𝑒 (1 + 𝑒 ) cos(3𝜔) + 𝐵1 − (1 + 3𝑒−2𝛼 + 5𝑒−4𝛼 + 7𝑒−6𝛼 ) sin(𝜔)+
⎪
⎪
⎪
)︂
⎪
⎪
⎪
⎪
⎪
+𝑒−2𝛼 (3 + 5𝑒−2𝛼 ) sin(3𝜔) = −(1 − 𝑒−2𝛼 )(3 + 4𝑒−2𝛼 + 3𝑒−4𝛼 ) cos(𝜔)+
⎪
⎪
⎨
+𝑒−2𝛼 (1 − 𝑒−2𝛼 ) cos(3𝜔)
⎪
⎪
⎪
(︂
)︂
⎪
⎪
⎪
−4𝛼
−4𝛼
⎪
⎪
𝐴1 4𝑒
− 2(1 − 3𝑒 ) cos(2𝜔) − 𝐵0 (1 + 4𝑒−2𝛼 + 𝑒−4𝛼 ) sin(2𝜔)+
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
+2𝐵1 (1 + 𝑒−2𝛼 )(1 + 3𝑒−2𝛼 ) sin(2𝜔) = −(1 − 𝑒−4𝛼 )(2 + cos(2𝜔))
⎪
⎪
⎪
⎪
⎪
⎪
−2𝛼
−2𝛼
−2𝛼
⎪
⎪
⎪𝐴1 (3 − 5𝑒 ) cos(𝜔) + 𝐵0 (1 + 𝑒 ) sin(𝜔) − 𝐵1 (3 + 5𝑒 ) sin(𝜔) =
⎪
⎪
⎪
⎪
⎩= −(1 − 𝑒−2𝛼 ) cos(𝜔)
Способом, рассмотренным в примере 1.1, находились спектральная плотность и
корреляционная функция у остальных примеров в данной работе.
10
1.2. Оценка коэффициентов стохастического уравнения
Рассмотрим 𝑋𝑡 — процесс АР(𝑛), удовлетворяющий стохастическому уравнению (1.2).
Считаем, что выбранный процесс стационарен, что означает, что все корни характери
стического полинома
𝑞(𝑦) =
𝑛
∑︁
𝑞𝑘 𝑦 𝑘 ,
𝑞0 = 1
(1.9)
𝑘=0
по модулю больше единицы.
Пусть 𝑋1 , 𝑋2 , . . . , 𝑋𝑁 — реализация рассматриваемого процесса. Вычислим выбо
рочные корреляции 𝑅𝑘𝑁 согласно формуле
𝑅𝑘𝑁
где 𝑥¯ = (1/𝑁 )
∑︀𝑁
𝑗=1
𝑁 −𝑘
1 ∑︁
=
(𝑋𝑗+𝑘 − 𝑥¯)(𝑋𝑗 − 𝑥¯),
𝑁 − 𝑘 𝑗=1
𝑘 = 0, 1, . . . , 𝑣,
(1.10)
𝑋𝑗 - выборочное среднее.
Далее будем оценивать параметры с помощью системы уравнений Юла — Уокера
⎧
⎪
⎪
𝑞0 𝑅1 + 𝑞1 𝑅0 + 𝑞2 𝑅1 + · · · + 𝑞𝑛 𝑅𝑛−1 = 0
⎪
⎪
⎪
⎪
⎪
⎪
⎨𝑞0 𝑅2 + 𝑞1 𝑅1 + 𝑞2 𝑅0 + · · · + 𝑞𝑛 𝑅𝑛−2 = 0
(1.11)
⎪
⎪
⎪. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
⎪
⎪
⎪
⎪
⎪
⎩𝑞 𝑅 + 𝑞 𝑅
0 𝑛
1 𝑛−1 + 𝑞2 𝑅𝑛−2 + · · · + 𝑞𝑛 𝑅0 = 0,
𝑝20 = 𝑅0 + 𝑞1 𝑅1 + · · · + 𝑞𝑛 𝑅𝑛 .
(1.12)
Учитывая, что 𝑞0 = 𝑞^0 = 1, можно решить данную систему относительно 𝑞1 , 𝑞2 , . . . , 𝑞𝑛 ,
подставив в неё известные выборочные корреляции и найти соответствующие оценки
𝑞^1 , 𝑞^2 , . . . , 𝑞^𝑛 . После этого воспользуемся формулой (1.12) и вычислим оценку для коэф
фициента 𝑝0 .
√︁
𝑝^0 = 𝑅0𝑁 + 𝑞^1 𝑅1𝑁 + · · · + 𝑞^𝑛 𝑅𝑛𝑁
Систему 1.11 можно продолжить
𝑞0 𝑅𝑛+𝑘 + 𝑞1 𝑅𝑛+𝑘−1 + · · · + 𝑞𝑛 𝑅𝑘 = 0,
𝑘 = 1, 2, . . . .
Это приводит к методу наименьших квадратов. Данный подход заключается в решении
системы линейных уравнений вида
𝑛
∑︁
𝑗=0
𝑞^𝑗
𝑣
∑︁
𝑘=1
𝑁
𝑁
𝑅𝑘−𝑡
= 0,
𝑅𝑘−𝑗
𝑞^0 = 1,
𝑡 = 1, . . . , 𝑛.
(1.13)
11
𝑁
𝑁
Заметим, что процесс 𝑋𝑡 стационарен, а значит, что 𝑅𝑚
= 𝑅−𝑚
, 1 ≤ 𝑚 ≤ 𝑣, и
для оценки коэффициентов 𝑞𝑗 необходимо располагать не менее, чем 𝑛 выборочными
𝑁
, то есть 𝑣 ≥ 𝑛. После нахождения решения этой системы мы получим
корреляциями 𝑅𝑚
оценки 𝑞^1 , 𝑞^2 , . . . , 𝑞^𝑛 .
В методе наименьших квадратов используется больше выборочных корреляций,
чем в системе Юла–Уокера. Это позволяет при построении модели учитывать больше
информации о процессе, что приводит к построению более эффективной модели. По
этой причине, в данной работе применялся метод наименьших квадратов.
Теорема 1.2 ([3]). Пусть 𝑋𝑡 — стационарный процесс АР(𝑛), удовлетворяющий сто
хастическому уравнению (1.2). Пусть с помощью уравнений Юла–Уокера вычислены
𝑞^1 , 𝑞^2 , . . . , 𝑞^𝑛 , а 𝑝^0 имеет вид
𝑝^0 = 𝜎
^
√︀
1 + 𝑞^1 𝜌^1 + · · · + 𝑞^𝑛 𝜌^𝑛 ,
𝜎^2 = 𝑅^0 ,
(1.14)
где коэффициенты корреляций 𝜌^𝑗 выражаются через корреляции 𝑅^𝑗 следующим обра
зом
𝜌^𝑗 =
𝑅^𝑗
,
𝑅^0
𝑗 = 1, 2 . . . , 𝑛.
(1.15)
Тогда при достаточно большом числе реализаций 𝑁
√
1
𝑁 (𝑞^𝑘 − 𝑞𝑘 ) ∼ 𝑁 (0, 𝑝20 2 𝑐𝑘𝑘 ),
𝜎
𝑘 = 1, 2 . . . , 𝑛,
(1.16)
где 𝜎 2 = 𝑅0 ,
⎛
⎜
⎜
⎜
𝑅=⎜
⎜
⎜
⎝
𝑅0
𝑅1
..
.
𝑅1
𝑅0
..
.
...
...
...
𝑅𝑛−1 𝑅𝑛−2 . . .
𝑅𝑛−1
⎞
⎛
1
𝜌1
...
⎟
⎜
⎟
⎜
⎜ 𝜌1
𝑅𝑛−2 ⎟
1
...
2
⎜
⎟
.. ⎟ = 𝜎 ⎜ ..
..
⎜ .
. ⎟
.
...
⎠
⎝
𝑅0
𝜌𝑘−1 𝜌𝑘−2 . . .
𝜌𝑘−1
⎞
⎟
⎟
𝜌𝑘−2 ⎟ −1
1
⎟
.. ⎟, 𝑅 = 𝜎 2
. ⎟
⎠
1
⎛
𝑐
...
⎜ 11
⎜ ..
⎜ . ...
⎝
𝑐𝑛1 . . .
⎞
𝑐1𝑛
⎟
.. ⎟
. ⎟.
⎠
𝑐𝑛𝑛
(1.17)
Следствие 1.1 ([3]). Если 𝑛 = 1, то 𝑞1 = −𝜌1 и
√︁
𝑋𝑡 − 𝜌1 𝑋𝑡−1 = 𝜎 1 − 𝜌21 𝜉𝑡 ,
√
𝑁 (𝜌^1 − 𝜌1 ) ∼ 𝑁 (0, (1 − 𝜌21 )).
(1.18)
В связи с этим, для проверки гипотезы
𝐻0 : 𝜌1 = 0
(1.19)
12
√
√
используется статистика 𝑇 = | 𝑁 𝜌^1 |: если | 𝑁 𝜌^1 | > 1.96, то с вероятностью 0.95 гипо
теза отвергается и считается, что 𝜌1 ̸= 0.
Следствие 1.2 ([3]). Если 𝑛 = 2, то
√︀
𝑋𝑡 + 𝑞1 𝑋𝑡−1 + 𝑞2 𝑋𝑡−2 = 𝜎 1 + 𝑞1 𝜌1 + 𝑞2 𝜌2 𝜉𝑡 ,
⎞
⎛
⎛
⎞
−𝜌1
1
1
𝜌
2
2
1
1
⎠,
𝑅−1 = 2 ⎝ 1−𝜌1 1−𝜌1 ⎠ ,
𝑅 = 𝜎2 ⎝
−𝜌
1
1
𝜎
𝜌1 1
1−𝜌21
1−𝜌21
)︂
(︂
√
1
,
𝑗 = 1, 2.
𝑁 (𝑞^𝑗 − 𝑞𝑗 ) ∼ 𝑁 0, (1 + 𝑞1 𝜌1 + 𝑞2 𝜌2 )
1 − 𝜌21
(1.20)
(1.21)
Согласно теореме 1.2 при заданном 𝑛 дисперсии параметров 𝑞^1 , 𝑞^2 , . . . , 𝑞^𝑛 имеют
порядок 𝑂( 𝑁1 ).
Таким образом, по наблюдениям над стационарным процессом мы смогли оценить
параметры стохастического уравнения, которому удовлетворяет рассматриваемый про
цесс. Соответственно, зная эти оценки, можно построить спектральную плотность и
корреляционную функцию для полученной модели.
Кроме того, построенную модель следует проверить на стационарность. Для этого,
согласно теореме 1.1, следует найти корни характеристического полинома у полученной
модели и убедиться, что они находятся вне единичного круга.
13
1.3. Определение порядка процесса авторегрессии
В предыдущем разделе был описан метод построения модели по ряду наблюдений,
а именно метод оценивания параметров. Но данная процедура зависит от порядка про
цесса авторегрессии 𝑛. В реальной жизни мы располагаем только лишь рядом наблю
дений, а про порядок процесса нам ничего неизвестно. В связи с этим для построения
эффективной модели возникает необходимость определения её порядка.
В работе [4] предлагается для решения этой проблемы воспользоваться частной
автокорреляционной функцией. Обозначим через 𝑞𝑘𝑗 𝑗 – тый коэффициент процесса
АР(𝑘). Последний коэффициент 𝑞𝑘𝑘 называется функцией частной автокорреляции.
Тогда систему уравнений Юла–Уокера 1.11 можно записать в виде
⎛
⎞ ⎛ ⎞
⎛ ⎞
1
𝜌1
𝜌2 . . . 𝜌𝑘−1
𝑞𝑘1
𝜌
⎜
⎟ ⎜ ⎟
⎜ 1⎟
⎜
⎟ ⎜ ⎟
⎜ ⎟
⎜ 𝜌1
⎜𝜌 ⎟
1
𝜌1 . . . 𝜌𝑘−2 ⎟ ⎜𝑞𝑘2 ⎟
⎜
⎟ · ⎜ ⎟ = − ⎜ 2⎟ ,
⎜ ..
⎜ .. ⎟
..
..
.. ⎟ ⎜ .. ⎟
⎜ .
⎜.⎟
.
.
...
. ⎟ ⎜ . ⎟
⎝
⎠ ⎝ ⎠
⎝ ⎠
𝜌𝑘−1 𝜌𝑘−2 𝜌𝑘−3 . . .
1
𝑞𝑘𝑘
𝜌𝑘
(1.22)
где коэффициенты корреляций 𝜌𝑗 вычисляются по формуле 1.15.
Решая эти уравнения последовательно при 𝑘 = 1, 2, 3, получаем
𝑞11 = −𝜌1 ,
⃒
⃒
⃒
⃒
⃒ 1 −𝜌1 ⃒
⃒
⃒
⃒
⃒
⃒𝜌1 −𝜌2 ⃒
𝜌 − 𝜌21
⃒ =− 2
,
𝑞22 = ⃒
⃒
⃒
1 − 𝜌21
⃒ 1 𝜌1 ⃒
⃒
⃒
⃒
⃒
⃒𝜌 1 1 ⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒ 1 𝜌1 −𝜌1 ⃒ ⧸︃ ⃒ 1 𝜌1 𝜌2 ⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒
𝑞33 = ⃒𝜌1 1 −𝜌2 ⃒
⃒𝜌1 1 𝜌1 ⃒ .
⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒
⃒𝜌2 𝜌1 −𝜌3 ⃒
⃒𝜌 2 𝜌 1 1 ⃒
(1.23)
В общем случае определитель в числителе состоит из тех же элементов, что и в
знаменателе, но последний столбец заменяется столбцом свободных членов из систе
мы 1.22.
Для процесса (1.2) порядка 𝑛 справедлива гипотеза
𝐻0 : 𝑞𝑛+1,𝑛+1 = 0.
(1.24)
Таким образом, определить порядок процесса можно с помощью частной автокор
реляционнной функции.
14
1.4. Выявление кратных корней характеристического полинома
В разделах 1.2 и 1.3 описаны методы построения эффективной модели по времен
ному ряду. Но при применении их на практике возникают некоторые сложности.
В случае, когда в знаменателе спектральной плотности нет кратных корней, по
частной автокорреляционной функции удается определить правильный порядок моде
ли. Затем, используя метод из раздела 1.2, находятся соответствующие оценки парамет
ров. Причем, построенный по этим оценкам полином 𝑞^(𝑧), будем иметь корни довольно
близкие к исходным.
В случае, когда знаменатель спектральной плотности обладает корнями с высоким
порядком кратности, этот подход не даёт положительного результата: по частной ав
токорреляционной функции порядок модели определить удается не всегда, и разница
между корнями полиномов 𝑞(𝑧) и 𝑞^(𝑧) оказывается существенной, что видно из при
мера, в котором рассматривается процесс 2.3. Поэтому возникла задача о нахождении
способа идентификации кратных корней.
В данной работе проводилось исследование поставленной задачи и был найден
метод, который позволяет выявить корни знаменателя спектральной плотности и опре
делить их кратность. Далее мы подробно опишем суть построенного метода.
Основная идея построенного метода заключается в последовательном исключении
составляющих компонент процесса.
Подход исключения используется в работе [5] для получения стационарного про
цесса. В. П. Носко утверждает, что если процесс авторегрессии 𝑋𝑡 нестационарен и
характеристический полином имеет единичный корень, то процесс вида 𝑌𝑡 = 𝑋𝑡 − 𝑋𝑡−1
может обладать свойством стационарности. Если стационарность не появилась, то для
процесса 𝑌𝑡 также необходимо выполнить исключение некоторой компоненты. Так, по
следовательно исключая мешающие компоненты, можно получить стационарный про
цесс.
Аналогичный подход применяется в разработанном алгоритме для выявления крат
ных корней знаменателя спектральной плотности.
Рассмотрим процесс АР(𝑛), причем 𝑛 = 2𝑘, 𝑘 > 1 и характеристический полином
знаменателя спектральной плотности 𝑞(𝑧) имеет пару комплексно–сопряжённых корней
15
кратности 𝑘. Представим эти корни в следующем виде
𝑧1,2 = 𝑒𝛼 (cos 𝜔 ± 𝑖 sin 𝜔).
Пусть имеется реализация этого процесса 𝑋1 , 𝑋2 , . . . , 𝑋𝑁 .
ШАГ 1. Реализуем процесс оценивания коэффициентов модели из раздела 1.2 по
наблюдениям 𝑋1 , 𝑋2 , . . . , 𝑋𝑁 , считая, что порядок модели 𝑛 = 2, и получим коэффици
(1)
(1)
енты 𝑞^1 , 𝑞^2 . Затем составим характеристический полином и найдем его корни.
(1)
(1)
𝑞^(1) (𝑧) = 1 + 𝑞^1 𝑧 + 𝑞^2 𝑧 2 ,
(1)
(1)
𝑧^1,2 = 𝑒𝛼 (cos 𝜔 (1) ± 𝑖 sin 𝜔 (1) ).
Сравним между собой 𝛼, 𝜔 и 𝛼(1) , 𝜔 (1) и увидим, что они будут достаточно близки
ми.
ШАГ 2. По исходным наблюдениям 𝑋1 , 𝑋2 , . . . , 𝑋𝑁 и вычисленным параметрам
(1)
(1)
𝑞^1 , 𝑞^2 составим процесс 𝑌𝑡 вида
(1)
(1)
𝑌𝑡 = 𝑋𝑡+2 + 𝑞^1 𝑋𝑡+1 + 𝑞^2 𝑋𝑡 .
Тогда вместо исходного ряда получаем новый 𝑌1 , . . . , 𝑌𝑁 −2 . Выполним ШАГ 1 с наблю
дениями 𝑌1 , . . . , 𝑌𝑁 −2 и получим
(2)
(2)
𝑞^(2) (𝑧) = 1 + 𝑞^1 𝑧 + 𝑞^2 𝑧 2 ,
(2)
(2)
𝑧^1,2 = 𝑒𝛼 (cos 𝜔 (2) ± 𝑖 sin 𝜔 (2) ).
И снова обнаружим, что 𝛼, 𝜔 и 𝛼(2) , 𝜔 (2) близки между собой.
ШАГ 3. Повторим последовательно ШАГ 2 еще (𝑘 − 2) раз, используя при этом
наблюдения и оценки коэффициентов, вычисленные на предыдущем шаге. Таким обра
зом, у нас образуются две последовательности пар
(1)
(1)
(2)
(2)
(𝑘)
(𝑘)
(^
𝑞1 , 𝑞^2 ), (^
𝑞1 , 𝑞^2 ), . . . , (^
𝑞1 , 𝑞^2 )
(𝛼(1) , 𝜔 (1) ), (𝛼(2) , 𝜔 (2) ), . . . , (𝛼(𝑘) , 𝜔 (𝑘) ).
При этом, все 𝛼(𝑗) будут близки к 𝛼, а все 𝜔 (𝑗) — к 𝜔.
ШАГ 4. Реализуем ШАГ 2 по наблюдениям и коэффициентам, полученным на
предшествующем шаге. Если получим два вещественных корня, у которых 𝜔 (𝑘+1) =
16
(𝑘+1)
0, а величины 𝛼1
(𝑘+1)
, 𝛼2
гораздо больше тех, что были на предыдущих шагах, то
необходимо вернуться и снова выполнить этот шаг, только при условии, что порядок
модели 𝑛 = 1. После такой замены произойдет выделение только одного вещественного
(𝑘+1)
корня вида ±𝑒𝛼
. Далее следует проверить гипотезу 1.19 с помощью неравенства
√ (𝑘+1)
| 𝑁 𝜌1
| < 1.96,
(𝑘+1)
𝜌1
= ±𝑒−𝛼
(𝑘+1)
Если неравенство выполняется, то это говорит о том, что с вероятностью 0.95 все основ
ные составляющие процесса мы уже выявили на предыдущих шагах, и остался только
так называемый белый шум. Иначе снова требуется провести исключение компонент
процесса, то есть повторение шага 2.
Итак, подведем итоги проделанных действий. Последовательно проделав проце
дуру формирования нового процесса и оценивания соответствующих параметров, мы
получили последовательность из 𝑘 пар вида (𝛼(1) , 𝜔 (1) ), (𝛼(2) , 𝜔 (2) ), . . . , (𝛼(𝑘) , 𝜔 (𝑘) ) такую,
что все первые компоненты близки к исходному значению 𝛼, все вторые компоненты —
к исходному значению 𝜔.
Известно, что при заданном 𝑛 оценивание параметров имеет порядок 𝑂( √1𝑁 ), где
𝑁 — число наблюдений процесса (см. [3]). Тогда если две пары (𝛼(𝑗) , 𝜔 (𝑗) ), (𝛼(𝑚) , 𝜔 (𝑚) ) по
обеим компонентам отличаются на величину порядка
√1 ,
𝑁
то можно считать, что пары
совпадают. Также согласно следствию 1.2, если для всех 𝑗 выполняются следующие
соотношения
𝑂(1)
|𝛼(𝑗) − 𝛼| < √ ,
𝑁
𝑂(1)
|𝜔 (𝑗) − 𝜔| < √
𝑁
то разницу между оценками и точными значениями можно считать несущественной.
Тогда можно сделать вывод, что на самом деле характеристический полином зна
менателя спектральной плотности имеет только одну пару комплексно–сопряженных
корней кратности 𝑘, что и требовалось выявить.
Заметим, что в случае, когда характеристический полином знаменателя спектраль
ной плотности 𝑞(𝑧) имеет вещественные корни, для выявления этих корней необходимо
проводить пошаговую процедуру с порядком модели 𝑛 = 1.
В следующем разделе приведена информация об оценке погрешностей параметров
𝛼
^и𝜔
^ у процесса АР(2).
17
1.5. Оценка погрешностей параметров 𝛼
^ и𝜔
^ у процесса АР(2)
Рассмотрим процесс АР(𝑛) в случае, когда 𝑛 = 2 и корни характеристического
полинома комплексно–сопряженные. Тогда в формуле 1.20
𝑞1 = −2𝑒−𝛼 cos(𝜔),
𝑞2 = 𝑒−2𝛼 .
По ряду наблюдений подбираем модель, в которой
𝑞^1 = −2𝑒−𝛼^ cos(^
𝜔 ),
𝛼
^ = 𝛼 + ℰ1 ,
𝑞^2 = 𝑒−2𝛼^ ,
𝜔
^ = 𝜔 + ℰ2 ,
где ℰ1 , ℰ2 — зависимые случайные величины с Eℰ1 = 0, Eℰ2 = 0. При большой выбор
ке дисперсии Dℰ1 = Eℰ12 и Dℰ2 = Eℰ22 малы. Для того, чтобы найти эти дисперсии,
нужно знать D𝑞^1 , D𝑞^2 , 𝑐𝑜𝑣(𝑞^1 , 𝑞^2 ). Покажем, что они удовлетворяют системе линейных
уравнений 1.25.
⎧
⎪
⎪
⎪
D𝑞^1 + 𝜌21 D𝑞^2 + 2𝜌1 𝑐𝑜𝑣(𝑞^1 , 𝑞^2 ) = D𝜌^1
⎪
⎪
⎨
𝜌21 D𝑞^1 + D𝑞^2 + 2𝜌1 𝑐𝑜𝑣(𝑞^1 , 𝑞^2 ) = D𝜌^2
⎪
⎪
⎪
⎪
⎪
⎩𝜌1 D𝑞^1 + 𝜌1 D𝑞^2 + (1 + 𝜌21 )𝑐𝑜𝑣(𝑞^1 , 𝑞^2 ) = 𝑐𝑜𝑣(𝜌^1 , 𝜌^2 )
1
𝑐𝑜𝑣(𝜌^1 , 𝜌^2 ) =
𝑁
(1.25)
(︂ 𝑁
−1
∑︁
(︂
)︂
(︂
)︂)︂
𝑁
−1
∑︁
|𝑘|
|𝑘|
𝜌𝑘 𝜌𝑘+1 1 −
𝜌𝑘−1 𝜌𝑘+2 1 −
+
𝑁
𝑁
𝑘=−𝑁 +1
𝑘=−𝑁 +1
Из уравнений Юла–Уокера для теоретических и выборочных корреляций можно
получить систему
⎧
⎪
⎨(𝑞^1 − 𝑞1 ) + 𝜌1 (𝑞^2 − 𝑞2 ) = −𝜌^1 + 𝜌1
(1.26)
⎪
⎩𝜌1 (𝑞^1 − 𝑞1 ) + (𝑞^2 − 𝑞2 ) = −𝜌^2 + 𝜌2
Вычисляя дисперсии левой и правой частей уравнений 1.26, получаются первые
два уравнения системы 1.25. Если перемножить левые части уравнений 1.26, взять от
них математическое ожидание и приравнять к математическому ожиданию от перемно
женных правых частей этих же уравнений, то получим третье уравнение системы 1.25.
Решение системы 1.25 будет иметь вид
D𝑞^1 =
𝐶1
,
𝑁
D𝑞^2 =
𝐶2
,
𝑁
𝑐𝑜𝑣(𝑞^1 , 𝑞^2 ) =
𝐶3
.
𝑁
18
Перейдем к вычислению оценок дисперсий величин ℰ1 , ℰ2 . Начнем с ℰ1 .
E𝑞^2 = E𝑒−2(𝛼+ℰ1 ) = 𝑒−2𝛼 E𝑒−2ℰ1 = 𝐴
D𝑞^2 = E𝑒−4(𝛼+ℰ1 ) − 𝐴2 = 𝑒−4𝛼 E𝑒−4ℰ1 − 𝐴2
(1.27)
Заметим, что
E𝑒−2ℰ1 = E(1 − 2ℰ1 + 2ℰ12 + . . . ) ≈ 1 + 2Eℰ12
𝐴 ≈ 𝑒−2𝛼 (1 + 2Eℰ12 ),
𝐴2 ≈ 𝑒−4𝛼 (1 + 4Eℰ12 ),
E𝑒−4ℰ1 = E(1 − 4ℰ1 + 8ℰ12 + . . . ) ≈ 1 + 8Eℰ12 .
Поэтому из формулы 1.27 следует, что
D𝑞^2 = 4𝑒−4𝛼 Eℰ12 =
Eℰ12 =
𝐶2
,
𝑁
𝐶2 4𝛼
𝑒 .
4𝑁
(1.28)
Теперь перейдем к ℰ2 .
E𝑞^1 = E(−2𝑒−(𝛼+ℰ1 ) cos(𝜔 + ℰ2 )) = 𝐵
𝐶1
(1.29)
𝑁
Воспользуемся формулами преобразования и разложения в ряд Тейлора тригоно
D𝑞^1 = E(4𝑒−2(𝛼+ℰ1 ) cos2 (𝜔 + ℰ2 )) − 𝐵 2 =
метрических функций.
𝐵 = −2𝑒−𝛼 E𝑒−ℰ1 (cos(𝜔) cos(ℰ2 ) − sin(𝜔) sin(ℰ2 )) =
)︂(︂
)︂
(︂
1 2
1 2
1 3
−𝛼
= −2𝑒 E 1 − ℰ1 + ℰ1 + . . .
cos(𝜔)(1 − ℰ2 + . . . ) − sin(𝜔)(ℰ2 − ℰ2 + . . . ) ≈
2
2
6
(︂
(︂
)︂
)︂
1
1
≈ −2𝑒−𝛼 cos(𝜔) 1 + Eℰ12 − Eℰ22 + sin(𝜔)Eℰ1 ℰ2
2
2
(︂
(︂
)︂
)︂
2
2
−2𝛼
2
2
𝐵 ≈ 4𝑒
cos (𝜔) 1 + Eℰ1 − Eℰ2 + sin(2𝜔)Eℰ1 ℰ2
E(4𝑒−2(𝛼+ℰ1 ) cos2 (𝜔 + ℰ2 )) = 2𝑒−2𝛼 E𝑒−2ℰ1 (1 + cos(2𝜔 + 2ℰ2 )) =
= 2𝑒−2𝛼 E𝑒−2ℰ1 (1 + cos(2𝜔) cos(2ℰ2 ) − sin(2𝜔) sin(2ℰ2 )) =
(︂
)︂(︂
−2𝛼
2
= 2𝑒 E 1 − 2ℰ1 + 2ℰ1 + . . .
1 + cos(2𝜔)(1 − 2ℰ22 + . . . ) −
)︂
8 3
− sin(2𝜔)(2ℰ2 − ℰ2 + . . . ) ≈
6
(︂
)︂
−2𝛼
2
2
2
2
≈ 4𝑒
cos (𝜔) + 2 cos (𝜔)Eℰ1 − cos(2𝜔)Eℰ2 + 2 sin(2𝜔)Eℰ1 ℰ2
19
В итоге получится
D𝑞^1 = 4𝑒−2𝛼 (cos2 (𝜔)Eℰ12 + sin2 (𝜔)Eℰ22 + sin(2𝜔)Eℰ1 ℰ2 ) =
𝐶1
.
𝑁
Далее, учитывая вычисленное значение Eℰ12 , составим систему из двух уравнений
следующего вида
⎧
⎪
⎨E𝑞^1 − 𝑞1 = 0
⎪
⎩D𝑞^1 =
𝐶1
𝑁
⇐⇒
⎧
⎪
⎨ 1 cos(𝜔)(Eℰ12 − Eℰ22 ) + sin(𝜔)Eℰ1 ℰ2 = 0
2
⎪
⎩cos2 (𝜔)Eℰ 2 + sin2 (𝜔)Eℰ 2 + sin(2𝜔)Eℰ1 ℰ2 =
1
2
𝐶1 2𝛼
𝑒
4𝑁
Решив построенную систему относительно Eℰ1 ℰ2 , Eℰ22 , получим
Eℰ22 =
Eℰ1 ℰ2 =
𝐶1 2𝛼
𝑒 ,
4𝑁
𝑒2𝛼 cos(𝜔)
(𝐶1 − 𝐶2 𝑒2𝛼 ).
8𝑁 sin(𝜔)
(1.30)
(1.31)
20
1.6. Моделирование процесса с дискретным параметром
Для исследования разработанного пошагового алгоритма требуется обладать на
блюдениями различных процессов. Для этого были изучены некоторые методы модели
рования случайных процессов с дробно–рациональными спектральными плотностями,
которые использовались для получения реализаций процессов. В данном разделе крат
ко представлен теоретический материал по данному вопросу, подробно изложенный
в [2].
Рассмотрим стационарный гауссовский процесс 𝑋𝑡 с E𝑋𝑡 = 0, удовлетворяющий
разностному уравнению (1.1). При этом будем считать, что полиномы
𝑝(𝑧) =
𝑚
∑︁
𝑝𝑘 𝑧 𝑘 ,
𝑘=0
𝑞(𝑧) =
𝑛
∑︁
𝑞𝑘 𝑧 𝑘
𝑘=0
имеют корни вне единичного круга.
Будем моделировать процесс 𝑋𝑡 при 𝑡 = 0, 1, 2, . . .
Для получения начального значения 𝑋0 вычислим вектор
𝑋 T = (𝑋−1 , . . . , 𝑋−𝑛 , 𝜉−1 , . . . , 𝜉−𝑚 )
и промоделируем независимую от компонент вектора 𝑋 гауссовскую случайную вели
чину 𝜉0 , а затем воспользуемся уравнением
𝑋0 = −𝑞1 𝑋−1 − · · · − 𝑞𝑛 𝑋−𝑛 + 𝑝0 𝜉0 + 𝑝1 𝜉−1 + · · · + 𝑝𝑚 𝜉−𝑚 .
Вектор 𝑋 имеет гауссовские компоненты
⎛
𝑅
𝑅1 . . .
⎜ 0
⎜
⎜ 𝑅1
𝑅0 . . .
⎜
⎜
⎜ ...
... ...
⎜
⎜
⎜𝑅𝑛−1 𝑅𝑛−2 . . .
𝐾𝑋 = ⎜
⎜
⎜ 𝑓0
0
...
⎜
⎜
⎜ 𝑓1
𝑓0 . . .
⎜
⎜
⎜ ...
... ...
⎝
𝑓𝑚−1 𝑓𝑚−2 . . .
(1.32)
и корреляционную матрицу 𝐾𝑋 = E(𝑋𝑋 T )
⎞
𝑅𝑛−1 𝑓0 𝑓1 . . . 𝑓𝑚−1
⎟
⎟
𝑅𝑛−2 0 𝑓0 . . . 𝑓𝑚−2 ⎟
⎟
⎟
... ... ... ... ... ⎟
⎟
⎟
𝑅0
0
0 ... ... ⎟
⎟.
⎟
0
1
0 ...
0 ⎟
⎟
⎟
0
0
1 ...
0 ⎟
⎟
⎟
... ... ... ... ... ⎟
⎠
...
0
0 ...
1
21
Взаимные корреляции 𝑓0 , 𝑓1 , . . . , 𝑓𝑚 и корреляции 𝑅0 , 𝑅1 , . . . , 𝑅𝑛 находятся из си
стем уравнений (1.33) и (1.34) соответственно.
⎧
⎪
⎪
𝑞0 𝑓 0 = 𝑝0
⎪
⎪
⎪
⎪
⎪
⎪
⎨𝑞 0 𝑓 1 + 𝑞 1 𝑓 0 = 𝑝 1
⎪
⎪
⎪
.....................
⎪
⎪
⎪
⎪
⎪
⎩𝑞 𝑓 + . . . + 𝑞 𝑓 = 𝑝
0 𝑚
𝑚 0
(1.33)
𝑚
⎧
⎪
⎪
𝑞0 𝑅0 + 𝑞1 𝑅1 + . . . + 𝑞𝑛 𝑅𝑛 = 𝑝0 𝑓0 + . . . + 𝑝𝑚 𝑓𝑚
⎪
⎪
⎪
⎪
⎪
⎪
⎪
𝑞0 𝑅1 + 𝑞1 𝑅0 + . . . + 𝑞𝑛 𝑅𝑛−1 = 𝑝1 𝑓0 + . . . + 𝑝𝑚 𝑓𝑚−1
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
....................................
⎪
⎪
⎨
𝑞0 𝑅𝑚 + 𝑞1 𝑅𝑚−1 + . . . + 𝑞𝑛 𝑅𝑛−𝑚 = 𝑝𝑚 𝑓0
⎪
⎪
⎪
⎪
⎪
⎪
𝑞0 𝑅𝑚+1 + 𝑞1 𝑅𝑚 + . . . + 𝑞𝑛 𝑅𝑛−𝑚−1 = 0
⎪
⎪
⎪
⎪
⎪
⎪
⎪
....................................
⎪
⎪
⎪
⎪
⎪
⎪
⎩𝑞 𝑅 + 𝑞 𝑅
0 𝑛
1 𝑛−1 + . . . + 𝑞𝑛 𝑅0 = 0
(1.34)
После вычисления элементов матрицы 𝐾𝑋 воспользуемся разложением Холецкого
для этой симметричной матрицы, а именно представим её в виде
𝐾𝑋 = 𝐿𝐿T ,
где 𝐿 — нижне–треугольная матрица.
Обозначим через 𝑊 (𝑘) независимые векторы
𝑊 (𝑘)T = (𝑤1 , 𝑤2 , . . . , 𝑤𝑘 ),
состоящие из 𝑘 независимых гауссовских стандартных компонент. Промоделировав 𝜉0
и 𝑊 (𝑚 + 𝑛), получим вектор 𝑋 = 𝐿𝑊 (𝑚 + 𝑛) и затем вычислим 𝑥0 по формуле (1.32).
Для получения каждого из последующих значений 𝑋𝑡 , 𝑡 = 1, 2, . . . , достаточно
промоделировать одно независимое от предыдущих 𝜉𝑡 и воспользоваться формулой
𝑋𝑡 = −𝑞1 𝑋𝑡−1 − · · · − 𝑞𝑛 𝑋𝑡−𝑛 + 𝑝0 𝜉𝑡 + 𝑝1 𝜉𝑡−1 + · · · + 𝑝𝑚 𝜉𝑡−𝑚 ,
(1.35)
поскольку к моменту времени 𝑡 остальные переменные, входящие в правую часть ра
венства (1.35), уже известны.
22
Глава 2
Практическое применение
2.1. Стационарный процесс с дискретным параметром
В ходе работы было проведено моделирование процесса с дискретным параметром.
В частности, сначала был рассмотрен 𝑋𝑡 — процесс АР(2), удовлетворяющий разност
ному уравнению
𝑋𝑡 − 1.553𝑋𝑡−1 + 0.854𝑋𝑡−2 = 0.287𝜉𝑡 .
(2.1)
Полином 𝑞(𝑧) имеет некратные корни 0.909 ± 0.587𝑖.
Спектральная плотность данного процесса имеет вид
Рис. 2.1. Спектральная плотность процесса 2.1
С помощью метода моделирования из раздела 1.6 была получена реализация дан
ного процесса, которая представлена на рисунке 2.14.
Затем рассматривался процесс 𝑋𝑡 чуть более сложный, а именно процесс АР(4),
удовлетворяющий разностному уравнению
𝑋𝑡 − 3.105𝑋𝑡−1 + 4.119𝑋𝑡−2 − 2.652𝑋𝑡−3 + 0.729𝑋𝑡−4 = 0.033𝜉𝑡 .
(2.2)
Полином 𝑞(𝑧) в данном случае имеет корни 0.909 ± 0.587𝑖, только кратности 2, а
спектральная плотность имеет вид
23
Рис. 2.2. Спектральная плотность процесса 2.2
С помощью метода моделирования из раздела 1.6 была получена реализация дан
ного процесса, которая представлена на рисунке 2.15.
Далее был рассмотрен еще более сложный процесс 𝑋𝑡 — процесс АР(6), удовлетво
ряющий разностному уравнению
𝑋𝑡 − 4.615𝑋𝑡−1 + 9.665𝑋𝑡−2 − 11.531𝑋𝑡−3 + 8.260𝑋𝑡−4 −
− 3.372𝑋𝑡−5 + 0.624𝑋𝑡−6 = 0.0028𝜉𝑡 .
(2.3)
В этом примере полином 𝑞(𝑧) имеет корни 0.9 ± 0.6𝑖 кратности 3, а спектральная
плотность изображена на рисунке 2.3.
Рис. 2.3. Спектральная плотность процесса 2.3
24
C помощью метода моделирования из раздела 1.6 была получена реализация этого
процесса. На рисунке 2.16 показан график изменения состояний процесса с течением
времени.
Заметим, что по графикам 2.14, 2.15, 2.16 уже видно, что процессы, обладающие
кратными корнями, отличаются от процессов с некратными корнями. На графиках 2.15
и 2.16 отображается периодичность, а на рисунке 2.14 она отсутствует.
Кроме того, в ходе работы были рассмотрены процессы с так называемым «меша
ющим воздействием». В качестве подобного воздействия, например, считается добав
ление к кратным корням некратного вещественного корня или пары некратных ком
плексно–сопряженных корней.
В частности, был рассмотрен процесс 𝑋𝑡 — процесс АР(6), удовлетворяющий раз
ностному уравнению
𝑋𝑡 − 4.808𝑋𝑡−1 + 10.167𝑋𝑡−2 − 12.041𝑋𝑡−3 + 8.404𝑋𝑡−4 −
− 3.278𝑋𝑡−5 + 0.559𝑋𝑡−6 = 0.0068𝜉𝑡 .
(2.4)
В этом примере полином 𝑞(𝑧) имеет корни 0.9 ± 0.6𝑖 кратности 2 и 1.13 ± 0.17𝑖
кратности 1, а спектральная плотность имеет вид
Рис. 2.4. Спектральная плотность процесса 2.4
C помощью метода моделирования из раздела 1.6 была получена реализация этого
процесса. На рисунке 2.17 показан график изменения состояний процесса с течением
времени.
25
Также был рассмотрен процесс 𝑋𝑡 — процесс АР(5), удовлетворяющий разностно
му уравнению
𝑋𝑡 − 3.986𝑋𝑡−1 + 6.873𝑋𝑡−2 − 6.336𝑋𝑡−3 + 3.121𝑋𝑡−4 − 0.664𝑋𝑡−5 = 0.016𝜉𝑡 .
(2.5)
В этом примере полином 𝑞(𝑧) имеет корни 0.9 ± 0.6𝑖 кратности 2 и 1.1 кратности
1, а спектральная плотность имеет вид
Рис. 2.5. Спектральная плотность процесса 2.5
C помощью метода моделирования из раздела 1.6 была получена реализация этого
процесса. На рисунке 2.18 показан график изменения состояний процесса с течением
времени.
Стоит отметить, что на графиках 2.17 и 2.18 периодичность отсутствует из–за
влияния «мешающего воздействия». По этой причине визуально определить наличие
кратных корней характеристического полинома не удается.
26
2.2. Подбор стационарной модели по методу Юла–Уокера
После того, как были получены реализации процессов, представленные в преды
дущем разделе, мы использовали их для оценки коэффициентов модели и нахождения
корней соответствующих характеристических полиномов. Сначала был применен метод
Юла–Уокера. Его идея заключается в том, что по ряду наблюдений необходимо вычис
лить частную автокорреляционную функцию и по ней определить порядок процесса
тем способом, который описан в разделе 1.3. Затем по вычисленному порядку процесса
строится соответствующая модель, параметры которой оцениваются способом из раз
дела 1.2. Далее в этом разделе будут описаны получившиеся результаты применения
этого метода.
Начнем с процесса АР(2), реализация которого представлена на рисунке 2.14. Эти
наблюдения были использованы для определения порядка процесса и получения оценок
коэффициентов модели. Поскольку данный пример очень простой: всего лишь второй
порядок процесса и нет кратных корней, поэтому достаточно легко был определен необ
ходимый порядок процесса, и после этого по ряду наблюдений была построена следую
щая модель
(2.1′ )
𝑋𝑡 − 1.574𝑋𝑡−1 + 0.875𝑋𝑡−2 = 0.269𝜉𝑡 .
Таблица 2.1. Результаты оценивания для процесса 2.1 по методу Юла–Уокера
Точные характеристики процесса 2.1
Полином 𝑞(𝑧)
1 − 1.553𝑧 + 0.854𝑧 2
Корни 𝑞(𝑧)
𝛼
𝑧1,2 = 0.909 ± 0.587𝑖 0.07887
𝜔
Период
0.57339
10.95796
𝜔
^
Период
0.57150,
10.9942
Характеристики1 процесса 2.1′
Полином 𝑞^(𝑧)
1 − 1.574𝑧 + 0.875𝑧 2
Корни 𝑞^(𝑧)
𝛼
^
𝑧^1,2 = 0.899 ± 0.578𝑖 0.06675,
𝜀1 = 0.01212
𝜀2 = 0.00189
Исходя из данных в таблице 2.1, можно сделать вывод, что используемый метод
достаточно хорошо оценил параметры процесса и выявил пару комплексно–сопряжен
ных корней близких к исходным. Чтобы удостоверится в правильности построенной
1
𝜀
(𝛼)
В этой таблице, а также во всех последующих, при представлении 𝛼
^, 𝜔
^ указываются величины
= |𝛼 − 𝛼
^ | и 𝜀(𝜔) = |𝜔 − 𝜔
^ |, обозначающие величины погрешности вычислений.
27
модели, на рисунке 2.6 представлены соответствующие коэффициенты корреляции2 .
Рис. 2.6. Корреляции и выборочные корреляции процесса 2.1. Корреляции процесса 2.1′
.
Далее оценивание проводилось для процесса АР(4), реализация которого представ
лена на рисунке 2.15. В данном примере в знаменателе спектральной плотности име
ются корни второй кратности. Тем не менее, был верно определен порядок процесса и
построена следующая модель
(2.2′ )
𝑋𝑡 − 2.983𝑋𝑡−1 + 3.827𝑋𝑡−2 − 2.385𝑋𝑡−3 + 0.641𝑋𝑡−4 = 0.0508𝜉𝑡 .
Таблица 2.2. Результаты оценивания для процесса 2.2 по методу Юла–Уокера
Точные характеристики процесса 2.2
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
0.07887
0.57339
10.95796
𝛼
^
𝜔
^
Период
1 − 2.983𝑧 + 3.827𝑧 2 − 𝑧^1,2 = 0.87 ± 0.567𝑖,
0.03697,
0.57645,
10.9,
2.385𝑧 3 + 0.641𝑧 4
𝜀1 = 0.0419, 𝜀1
1 − 3.105𝑧 + 4.119𝑧 2 − 𝑧1,2 = 0.909 ± 0.587𝑖
2.652𝑧 3 + 0.729𝑧 4
кратности 2
Характеристики процесса 2.2′
Полином 𝑞^(𝑧)
Корни 𝑞^(𝑧)
(𝛼)
𝑧^3,4 = 0.99 ± 0.685𝑖
2
(𝜔)
= 0.0031,
0.18530,
0.60510,
(𝛼)
𝜀2
(𝜔)
𝜀2
= 0.1064
10.38
= 0.0317
В данной работе моделируются процессы с единичной дисперсией, поэтому корреляции совпа
дают с коэффициентами корреляции.
28
Результат вычислений для данной модели представлен в таблице 2.2. По этим ре
зультатам можно отметить следующее: рассматриваемый метод позволил правильно
оценить коэффициенты модели и определить корни характеристического полинома.
Аналогичные выводы можно сделать из рисунка 2.7, на котором представлены как
корреляции исходного процесса, так и корреляции построенной модели.
Рис. 2.7. Корреляции и выборочные корреляции процесса 2.2. Корреляции процесса 2.2′
Теперь перейдем к рассмотрению еще более сложного процесса, а именно к про
цессу АР(6), реализация которого представлена на рисунке 2.16. В этом примере зна
менатель спектральной плотности обладает корнями третьей кратности. Наблюдения,
изображенные на рисунке 2.16, были использованы для определения порядка процесса и
построения соответствующей модели. В данном примере метод Юла–Уокера установил,
что порядок процесса равен 7, и была построена модель АР(7), которая удовлетворяет
стохастическому уравнению
𝑋𝑡 − 3.315𝑋𝑡−1 + 4.514𝑋𝑡−2 − 2.291𝑋𝑡−3 − 1.045𝑋𝑡−4 +
+ 2.182𝑋𝑡−5 − 1.241𝑋𝑡−6 + 0.296𝑋𝑡−7 = 0.016𝜉𝑡 .
(2.3′ )
В результате вычислений получилось, что рассматриваемый метод не выявил все
кратные корни: обнаружились две пары корней близкие к исходной и одна пара, сильно
отличающаяся от них. Поскольку две пары корней удалось выявить, мы попробовали
увеличить порядок модели с целью обнаружить еще одну пару. Поэтому была построена
29
модель АР(8), которая удовлетворяют стохастическому уравнению
𝑋𝑡 − 3.0.6𝑋𝑡−1 + 3.519𝑋𝑡−2 − 0.927𝑋𝑡−3 − 1.105𝑋𝑡−4 +
+ 0.116𝑋𝑡−5 + 1.449𝑋𝑡−6 − 1.237𝑋𝑡−7 + 0.359𝑋𝑡−8 = 0.016𝜉𝑡 .
(2.3′′ )
Таблица 2.3. Результаты оценивания для процесса 2.3 по методу Юла–Уокера
Точные характеристики процесса 2.3
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 4.615𝑧 + 9.665𝑧 2 −
𝑧1,2 = 0.9 ± 0.6𝑖
0.07706
0.58871
10.6728
11.531𝑧 3 + 8.26𝑧 4 −
кратности 3
3.372𝑧 5 + 0.624𝑧 6
Характеристики процесса 2.3′
Полином 𝑞^(𝑧)
Корни 𝑞^(𝑧)
𝛼
^
𝜔
^
Период
1 − 3.315𝑧 + 4.514𝑧 2 −
𝑧^1,2 = 0.897 ± 0.542𝑖,
0.04715,
0.54354,
11.56,
(𝛼)
2.291𝑧 3 − 1.045𝑧 4 +
2.182𝑧 5 − 1.241𝑧 6 +
𝜀1
𝑧^3,4 = 0.836 ± 0.6𝑖,
0.02877,
(𝛼)
0.296𝑧 7
𝜀2
𝑧^5,6 = 1.001 ± 1.123𝑖,
= 0.05,
0.4087,
(𝛼)
𝜀3
𝑧^7 = −1.28
= 0.03,
= 0.3,
0.24707,
(𝛼)
𝜀4
= 0.17
(𝜔)
𝜀1
= 0.05,
0.0.62278,
(𝜔)
𝜀2
= 0.03,
0.84275,
(𝜔)
𝜀3
(𝜔)
7.46,
= 0.25,
3.14159,
𝜀4
10.09,
2
= 2.6
Характеристики процесса 2.3′′
Полином 𝑞^(𝑧)
Корни 𝑞^(𝑧)
𝛼
^
𝜔
^
Период
1 − 3.06𝑧 + 3.519𝑧 2 −
𝑧^1,2 = 0.842 ± 0.602𝑖,
0.03463,
0.62074,
10.12,
(𝛼)
0.927𝑧 3 − 1.105𝑧 4 +
0.116𝑧 5 + 1.448𝑧 6 −
𝜀1
𝑧^3,4 = 0.906 ± 0.55𝑖,
0.05822,
(𝛼)
1.237𝑧 7 + 0.359𝑧 8
𝜀2
𝑧^5,6 = 0.994 ± 0.763𝑖,
= 0.02,
0.22587,
(𝛼)
𝜀3
𝑧^7,8 = −1.022 ± 0.653
= 0.04,
= 0.15,
(𝜔)
𝜀1
= 0.03,
0.0.54563,
(𝜔)
𝜀2
= 0.04,
0.65446,
(𝜔)
𝜀3
9.6,
= 0.07,
0.19282,
2.57338,
(𝛼)
𝜀4
(𝜔)
𝜀4
= 0.12
11.52,
2.44
= 1.98
Результаты проделанных вычислений для двух моделей представлены в табли
це 2.3. По этим результатам можно отметить следующее: при построении модели с
30
более высоким порядком удалось выявить три пары корней, которые близки к исход
ной, но также в этой модели присутствует одна пара корней, которая отличается от
остальных. Кроме того, обе построенные модели имеют порядок выше, чем у исходно
го процесса. В данном случае правильно определить порядок модели и оценить корни
характеристического полинома знаменателя спектральной плотности не удалось.
Тем не менее, корреляции построенной модели 2.3′′ оказываются довольно близки
ми к корреляциям исходного процесса. Подобный вывод можно сделать, проанализиро
вав рисунок 2.8.
Рис. 2.8. Корреляции и выборочные корреляции процесса 2.3. Корреляции процесса 2.3′′
Также в ходе работы были построены модели для процессов с «мешающим воздей
ствием».
В частности, рассмотрим процесс АР(6), реализация которого представлена на
рисунке 2.17. В этом примере в знаменателе спектральной плотности присутствует
некратные корни, выступающие в качестве «мешающего воздействия». Наблюдения,
изображенные на рисунке 2.17, были использованы для определения порядка процесса
и построения соответствующей модели АР(6), которая удовлетворяют стохастическому
уравнению
𝑋𝑡 − 2.714𝑋𝑡−1 + 2.649𝑋𝑡−2 − 0.627𝑋𝑡−3 − 0.575𝑋𝑡−4 +
+ 0.333𝑋𝑡−5 − 0.0029𝑋𝑡−6 = 0.072𝜉𝑡 .
(2.4′ )
Результат вычислений для данной модели представлен в таблице 2.4. По этим ре
зультатам можно отметить следующее: рассматриваемый метод не выявил все кратные
корни: обнаружились две пары корней близкие к исходным и два вещественных кор
31
ня. В данном случае также не удалось правильно оценить корни характеристического
полинома знаменателя спектральной плотности.
Таблица 2.4. Результаты оценивания для процесса 2.4 по методу Юла–Уокера
Точные характеристики процесса 2.4
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 4.808𝑧 + 10.167𝑧 2 −
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.09372,
0.58497,
10.74,
12.041𝑧 3 + 8.404𝑧 4 −
ности 2,
3.278𝑧 5 + 0.559𝑧 6
𝑧3,4 = 1.13 ± 0.17𝑖
0.13289
0.14639
42.92
кратности 1
Характеристики процесса 2.4′
Полином 𝑞^(𝑧)
Корни 𝑞^(𝑧)
𝛼
^
𝜔
^
Период
1 − 2.713𝑧 + 2.649𝑧 2 −
𝑧^1,2 = 0.863 ± 0.567𝑖,
0.0316,
0.58135,
10.81,
(𝛼)
0.627𝑧 3 − 0.575𝑧 4 +
0.333𝑧 5 − 0.029𝑧 6
𝜀1
𝑧^3,4 = 1.309 ± 0.236𝑖,
𝑧^5,6 = −2.01, 9
= 0.06,
(𝜔)
𝜀1
= 0.004,
0.28487,
0.17824,
(𝛼)
𝜀2
(𝜔)
𝜀2
= 0.15,
= 0.03,
0.69847,
3.14159,
(𝛼)
𝜀3
(𝜔)
𝜀3
= 0.56
35.25,
2
= 2.56
Тем не менее, корреляции построенной модели оказываются довольно близкими
к корреляциям исходного процесса. Подобный вывод можно сделать, проанализировав
рисунок 2.9.
Рис. 2.9. Корреляции и выборочные корреляции процесса 2.4. Корреляции процесса 2.4′
32
Также рассмотрим процесс АР(5), реализация которого представлена на рисун
ке 2.18. В этом примере в знаменателе спектральной плотности присутствует в каче
стве «мешающего воздействия» некратный вещественный корень. Наблюдения, изобра
женные на рисунке 2.18, были использованы для построения соответствующей модели
АР(5), которая удовлетворяют стохастическому уравнению
𝑋𝑡 − 3.763𝑋𝑡−1 + 6.109𝑋𝑡−2 − 5.256𝑋𝑡−3 + 2.392𝑋𝑡−4 − 0.463𝑋𝑡−5 = 0.018𝜉𝑡 .
(2.5′ )
Таблица 2.5. Результаты оценивания для процесса 2.5 по методу Юла–Уокера
Точные характеристики процесса 2.5
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 3.986𝑧 + 6.873𝑧 2 −
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.08137,
0.59155,
10.62156
6.336𝑧 3 + 3.121𝑧 4 −
ности 2, 𝑧3 = 1.1
0.09548
0
0.664𝑧 5
кратности 1
Характеристики процесса 2.5′
Полином 𝑞^(𝑧)
Корни 𝑞^(𝑧)
𝛼
^
𝜔
^
Период
1 − 3.763𝑧 + 6.109𝑧 2 −
𝑧^1,2 = 0.855 ± 0.575𝑖,
0.02959,
0.59186,
10.62,
(𝛼)
5.256𝑧 3 + 2.392𝑧 4 −
0.463𝑧 5
𝜀1
𝑧^3,4 = 1.093 ± 0.637𝑖,
0.23482,
(𝛼)
𝜀2
𝑧^5 = 1.273
= 0.05,
= 0.15,
0.24151,
(𝛼)
𝜀3
(𝜔)
𝜀1
= 0.0003,
0.52751,
(𝜔)
𝜀2
= 0.060,
(𝜔)
0, 𝜀3
11.91
=0
= 0.15
Результат вычислений для данной модели представлен в таблице 2.5. По этим ре
зультатам можно отметить следующее: рассматриваемый метод выявил «мешающий»
вещественный корень и две пары комплексно–сопряженных корней, близких к исходной.
Также корреляции построенной модели оказываются близкими к корреляциям исход
ного процесса, такой вывод следует из рисунка 2.10.
33
Рис. 2.10. Корреляции и выборочные корреляции процесса 2.5. Корреляции процесса 2.5′
34
2.3. Подбор стационарной модели с помощью пошагового метода
В предыдущем разделе были представлены результаты построения моделей по ме
тоду Юла–Уокера. И на примере процесса АР(6) мы убедились в том, что этот метод
не гарантирует как точного оценивания параметров процесса, так и идентификации
кратных корней характеристического полинома.
По этой причине к рассматриваемым примерам был применен пошаговый алго
ритм выявления кратных корней, описанный в разделе 1.4. Далее будут показаны по
лучившиеся результаты.
Начнем с процесса АР(2), реализация которого представлена на рисунке 2.14. Для
данного процесса пошаговый алгоритм совпадает с методом Юла–Уокера. Поэтому ре
зультаты оценивания получились идентичными. Но для полного убеждения, что по
шаговый алгоритм работает правильно, мы проделали дополнительный шаг (ШАГ 2).
Результат представлен в следующей таблице.
Таблица 2.6. Результаты оценивания для процесса 2.1 по пошаговому алгоритму
Точные характеристики процесса 2.1
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
1 − 1.553𝑧 + 0.854𝑧 2
𝑧1,2 = 0.909 ± 0.587𝑖
𝛼
𝜔
Период
0.07887
0.57339
10.95796
Результат оценивания при 𝑛 = 2, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
1 − 1.574𝑧 + 0.875𝑧 2
𝑧1,2 = 0.899 ± 0.578𝑖
(1)
𝛼(1)
𝜔 (1)
Период
0.06675,
0.57150,
10.9942
𝜀1 = 0.01212
𝜀2 = 0.00189
Результат оценивания при 𝑛 = 1, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
1 + 0.0094𝑧 − 0.0067𝑧 2
𝑧1,2 = 12.9
(2)
𝛼(2)
𝜔 (2)
2.55766
0
Период
—
По данной таблице видно, что исследуемый алгоритм привел к нужному результа
ту: удалось выявить два комплексно–сопряженных корня, близких к исходным, а после
ее исключения остался только белый шум.
Далее оценивание по пошаговому алгоритму проводилось для процесса АР(4), ре
ализация которого представлена на рисунке 2.15. В данном примере в знаменателе
спектральной плотности имеются корни второй кратности. Поэтому при использовании
35
описанного в разделе 1.4 алгоритма потребовалось проделать немного больше действий,
чем для первого рассматриваемого примера. По ряду наблюдений 2.15 был проведен
пошаговый алгоритм, результаты которого представлены в таблице 2.7.
Таблица 2.7. Результаты оценивания для процесса 2.2 по пошаговому алгоритму
Точные характеристики процесса 2.2
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
1 − 3.105𝑧 + 4.119𝑧 2 − 𝑧1,2 = 0.909 ± 0.587𝑖
2.652𝑧 3 + 0.729𝑧 4
𝛼
𝜔
Период
0.07887
0.57339
10.95796
кратности 2
Результат оценивания при 𝑛 = 2, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
1 − 1.637𝑧 + 0.9495𝑧 2
𝑧1,2 = 0.862 ± 0.557𝑖
(1)
𝛼(1)
𝜔 (1)
0.02589,
0.57380,
(𝛼)
𝜀1 = 0.05
Период
10.95
(𝜔)
𝜀1 = 0.0004
Результат оценивания при 𝑛 = 2, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
1 − 1.424𝑧 + 0.721𝑧 2
𝑧1,2 = 0.988 ± 0.642𝑖
(2)
𝛼(2)
𝜔 (2)
0.16383,
0.57590,
(𝛼)
𝜀2 = 0.08
Период
10.91
(𝜔)
𝜀2 = 0.003
Результат оценивания при 𝑛 = 1, ШАГ 3
Полином 𝑞^(3) (𝑧)
Корни 𝑞^(3) (𝑧)
1 − 0.042𝑧
𝑧1 = 23.599
(3)
𝛼(3)
𝜔 (3)
3.16122
0
Период
—
Рис. 2.11. Корреляции и выборочные корреляции процесса 2.2. Корреляции построенной мо
дели
36
По полученным результатам можно отметить следующее: описанный алгоритм
смог выявить две пары комплексно–сопряженных корней, которые близки между собой
и близки к исходной паре корней. Аналогичный вывод можно сделать из рисунка 2.11.
Теперь перейдем к рассмотрению еще более сложного процесса, а именно к про
цессу АР(6), реализация которого представлена на рисунке 2.16. В этом примере зна
менатель спектральной плотности обладает корнями третьей кратности. С помощью
пошагового алгоритма для данного процесса были получены следующие вычисления.
Таблица 2.8. Результаты оценивания для процесса 2.3 по пошаговому алгоритму
Точные характеристики процесса 2.3
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1−4.615𝑧+9.665𝑧 2 −11.531𝑧 3 +
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.07706
0.58871
10.6728
8.26𝑧 4 − 3.372𝑧 5 + 0.624𝑧 6
ности 3
𝜔 (1)
Период
0.58724,
10.6995
Результат оценивания при 𝑛 = 2, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
𝛼(1)
1 − 1.641𝑧 + 0.972𝑧 2
𝑧1,2 = 0.844 ± 0.562𝑖 0.01434,
(1)
(𝛼)
𝜀1 = 0.06
(𝜔)
𝜀1 = 0.001
Результат оценивания при 𝑛 = 2, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
𝛼(2)
1 − 1.563𝑧 + 0.862𝑧 2
𝑧1,2 = 0.907 ± 0.582𝑖 0.07433,
(2)
(𝛼)
𝜀2 = 0.003
𝜔 (2)
Период
0.57043,
11.0148
(𝜔)
𝜀2 = 0.02
Результат оценивания при 𝑛 = 2, ШАГ 3
Полином 𝑞^(3) (𝑧)
Корни 𝑞^(3) (𝑧)
𝛼(3)
1 − 1.343𝑧 + 0.672𝑧 2
𝑧1,2 = 0.999 ± 0.699𝑖 0.19874,
(3)
(𝛼)
𝜀3 = 0.12
𝜔 (3)
Период
0.6109,
10.2851
(𝜔)
𝜀3 = 0.02
Результат оценивания при 𝑛 = 1, ШАГ 4
Полином 𝑞^(4) (𝑧)
Корни 𝑞^(4) (𝑧)
1 − 0.046𝑧
𝑧1 = 21.74087
(4)
𝛼(4)
𝜔 (4)
3.07919
0
Период
—
Для рассматриваемого процесса АР(6) по таблице 2.8 можно сделать следующие
выводы:
∙ Пошаговый алгоритм оценивания дает хорошие результаты: алгоритм обнаружил
37
три пары комплексно–сопряженных корней, которые действительно близки к ис
ходной паре;
∙ Оценивание коэффициентов модели по методу Юла–Уокера не позволяет выявить
кратные корни: обнаруживаются две пары корней близкие к исходной и одна пара,
сильно отличающаяся от них. В данном случае верно оценить корни характери
стического полинома знаменателя спектральной плотности не удалось.
Также обратим внимание на сходство корреляционных функций процесса 2.3 и
построенной модели, изображенных на рисунке 2.12.
Рис. 2.12. Корреляции и выборочные корреляции процесса 2.3. Корреляции построенной мо
дели
Рассмотрим процесс АР(6), реализация которого представлена на рисунке 2.17. В
этом примере знаменатель спектральной плотности обладает корнями первой и второй
кратностей. С помощью пошагового алгоритма для данного процесса были проведены
вычисления, результаты которых приведены в таблице 2.9.
По полученным результатам можно отметить следующее: описанный алгоритм
смог выявить две пары комплексно–сопряженных корней, которые близки к исходной
паре корней кратности 2, и одну пару комплексно–сопряженных корней, которая вы
ступала в роли «мешающего воздействия».
Также обратим внимание на сходство корреляционных функций процесса 2.4 и по
лученной модели, изображенных на рисунке 2.13.
38
Таблица 2.9. Результаты оценивания для процесса 2.4 по пошаговому алгоритму
Точные характеристики процесса 2.4
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 4.808𝑧 + 10.167𝑧 2 −
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.09372,
0.58497,
10.74,
12.041𝑧 3 + 8.404𝑧 4 − 3.278𝑧 5 +
ности 2, 𝑧3,4 = 1.13 ±
0.13289
0.14639
42.92
0.559𝑧 6
0.17𝑖 кратности 1
𝜔 (1)
Период
0.53511,
11.74
Результат оценивания при 𝑛 = 2, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
𝛼(1)
1 − 1.606𝑧 + 0.871𝑧 2
𝑧1,2 = 0.922 ± 0.546𝑖 0.06909,
(1)
(𝛼)
𝜀1 = 0.02
(𝜔)
𝜀1 = 0.05
Результат оценивания при 𝑛 = 2, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
1 − 1.544𝑧 + 0.648𝑧 2
𝑧1,2 = 1.19 ± 0.35𝑖
(2)
𝛼(2)
𝜔 (2)
Период
0.21703,
0.28758,
21.85
(𝛼)
𝜀2 = 0.08
(𝜔)
𝜀2 = 0.14
Результат оценивания при 𝑛 = 2, ШАГ 3
Полином 𝑞^(3) (𝑧)
Корни 𝑞^(3) (𝑧)
1 − 1.499𝑧 + 0.822𝑧 2
𝑧1,2 = 0.912 ± 0.62𝑖
(3)
𝛼(3)
𝜔 (3)
Период
0.09794,
0.5974,
10.52
(𝛼)
𝜀3 = 0.004
(𝜔)
𝜀3 = 0.01
Результат оценивания при 𝑛 = 1, ШАГ 4
Полином 𝑞^(4) (𝑧)
Корни 𝑞^(4) (𝑧)
1 − 0.297𝑧
𝑧1 = 3.364
(4)
𝛼(4)
𝜔 (4)
Период
1.21313
0
—
Рис. 2.13. Корреляции и выборочные корреляции процесса 2.4. Корреляции построенной мо
дели.
39
Рассмотрим процесс АР(5), реализация которого представлена на рисунке 2.16. В
этом примере знаменатель спектральной плотности обладает корнями первой и второй
кратностей.
Для данного процесса был реализован пошаговый алгоритм в двух вариантах: в
первом — сначала выявление вещественного корня, а затем выявление кратных корней,
во втором — наоборот. Результаты вычислений представлены в таблицах 2.10 и 2.11.
Таблица 2.10. Результаты оценивания для процесса 2.5 по пошаговому алгоритму (вариант 1)
Точные характеристики процесса 2.5
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 3.986𝑧 + 6.873𝑧 2 − 6.336𝑧 3 +
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.08137,
0.59155,
10.62156
3.121𝑧 4 − 0.664𝑧 5
ности 2, 𝑧3 = 1.1
0.09548
0
кратности 1
Результат оценивания при 𝑛 = 1, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
1 − 0.765𝑧
𝑧1 = 1.307
(1)
𝛼(1)
𝜔 (1)
0.26786,
0,
(𝛼)
𝜀1
(𝜔)
𝜀1
= 0.17
Период
—
=0
Результат оценивания при 𝑛 = 2, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
𝛼(2)
𝜔 (2)
1 − 1.626𝑧 + 0.949𝑧 2
𝑧1,2 = 0.856 ± 0.565𝑖 0.02595,
(2)
(𝛼)
𝜀2
= 0.06
Период
0.58355,
(𝜔)
𝜀2
10.77
= 0.008
Результат оценивания при 𝑛 = 2, ШАГ 3
Полином 𝑞^(3) (𝑧)
Корни 𝑞^(3) (𝑧)
1 − 1.4004𝑧 + 0.649𝑧 2
𝑧1,2 = 1.08 ± 0.61𝑖
(3)
𝛼(3)
𝜔 (3)
0.21602,
0.51754,
(𝛼)
𝜀3
(𝜔)
𝜀3
= 0.13
Период
12.14
= 0.07
Результат оценивания при 𝑛 = 1, ШАГ 4
Полином 𝑞^(4) (𝑧)
Корни 𝑞^(4) (𝑧)
1 − 0.2497𝑧
𝑧1 = 4.005
(4)
𝛼(4)
𝜔 (4)
1.38751
0
Период
—
40
Таблица 2.11. Результаты оценивания для процесса 2.5 по пошаговому алгоритму (вариант 2)
Точные характеристики процесса 2.5
Полином 𝑞(𝑧)
Корни 𝑞(𝑧)
𝛼
𝜔
Период
1 − 3.986𝑧 + 6.873𝑧 2 − 6.336𝑧 3 +
𝑧1,2 = 0.9 ± 0.6𝑖 крат
0.08137,
0.59155,
10.62156
3.121𝑧 4 − 0.664𝑧 5
ности 2, 𝑧3 = 1.1
0.09548
0
кратности 1
Результат оценивания при 𝑛 = 2, ШАГ 1
Полином 𝑞^(1) (𝑧)
Корни 𝑞^(1) (𝑧)
1 − 1.521𝑧 + 0.616𝑧 2
𝑧1,2 = 0.903 ± 0.55𝑖
(1)
𝛼(1)
𝜔 (1)
0.0.05611,
0.54699,
(𝛼)
𝜀1 = 0.03
Период
11.49
(𝜔)
𝜀1 = 0.04
Результат оценивания при 𝑛 = 2, ШАГ 2
Полином 𝑞^(2) (𝑧)
Корни 𝑞^(2) (𝑧)
1 − 1.626𝑧 + 0.949𝑧 2
𝑧1,2 = 1.24 ± 0.0.31𝑖
(2)
𝛼(2)
𝜔 (2)
0.24248,
0.24687,
(𝛼)
𝜀2 = 0.16
Период
25.45
(𝜔)
𝜀2 = 0.34
Результат оценивания при 𝑛 = 1, ШАГ 3
Полином 𝑞^(3) (𝑧)
Корни 𝑞^(3) (𝑧)
1 − 0.632𝑧
𝑧1 = 1.581
(3)
𝛼(3)
𝜔 (3)
0.45822,
0,
(𝛼)
𝜀3 = 0.36
Период
—
(𝜔)
𝜀3 = 0
Результат оценивания при 𝑛 = 1, ШАГ 4
Полином 𝑞^(4) (𝑧)
Корни 𝑞^(4) (𝑧)
1 − 0.297𝑧
𝑧1 = 3.364
(4)
𝛼(4)
𝜔 (4)
0.70236
0
Период
—
По приведенным таблицам можно сделать следующие выводы:
∙ В первом варианте удалось правильно определить две пары комплексно–сопря
женных корней, которые оказались достаточно близки к исходной кратной паре,
и выявить «мешающий» вещественный корень;
∙ Во втором варианте найденные пары комплексно–сопряженных корней как несхо
жи между собой, так и отличаются от исходной пары корней второй кратности,
разница между «мешающим» вещественным корнем и полученной оценкой суще
ственна.
41
Таким образом, можно сказать, что для построения эффективной модели в по
добных ситуациях необходимо сначала выявить и исключить из рассмотрения «меша
ющий» вещественный корень, а затем оценивать и идентифицировать кратные ком
плексно–сопряженные корни.
42
2.4. Сравнительная характеристика методов
В ходе работы на довольно разнообразных примерах были исследованы два спосо
ба оценивания параметров стационарных процессов: метод Юла–Уокера и пошаговый
алгоритм исключения. По полученным результатам можно провести сравнительную
характеристику данных методов.
Метод Юла–Уокера достаточно точно определяет порядок процесса и строит со
ответствующую модель в случаях, когда знаменатель спектральной плотности имеет
корни кратности не выше второй. Такому выводу соответствуют результаты вычисле
ний, приведенных в таблицах 2.1, 2.2. В случаях третьей кратности и выше, а также
в случаях с «мешающим» воздействием правильное оценивание корней провести не
удается, это показывают результаты из таблиц 2.3 и 2.4.
Пошаговый алгоритм исключения позволяет правильно оценить параметры про
цесса по временному ряду для построения эффективной модели. Кроме того, алгоритм
дает возможность идентификации количества кратных корней знаменателя спектраль
ной плотности и определения их кратности. Также найденный алгоритм применим в
случаях с высоким порядком кратности корней и при наличии «мешающего» воздей
ствия. Подобный вывод можно сделать на основе проведенных вычислений для рас
сматриваемых процессов, представленных в таблицах 2.7, 2.8, 2.9 и 2.10.
Таблица 2.12. Величины погрешности вычислений
Процесс
Метод Юла–Уокера
Пошаговый алгоритм
𝜀(𝛼)
𝜀(𝜔)
𝜀(𝛼)
𝜀(𝜔)
2.2, кратность 2
0.04, 0.11
0.003, 0.03
0.05, 0.08
0.0004, 0.003
2.3, кратность 3
0.03, 0.05, 0.3 0.05, 0.03, 0.25 0.06, 0.003, 0.12
2.4, кратность 2 и 1 0.06, 0.15,
(компл. + компл.)
0.56
2.5, кратность 2 и 1 0.05, 0.15,
(компл.+ веществ.)
0.15
0.001, 0.02, 0.02
0.004, 0.03,
0.02, 0.08,
0.05, 0.14,
2.56
0.004
0.01
0.0003, 0.06,
0.06, 0.13,
0.008, 0.07,
0
0.17
0
На рассмотренных примерах видно, что пошаговый алгоритм оценивает лучше,
чем метод Юла–Уокера: при пошаговом алгоритме полученные оценки являются бо
лее точными и более близкими к истинным значениям параметров. То есть величины
43
погрешности вычислений меньше в случае пошагового оценивания, чем в случае оцени
вания по методу Юла–Уокера. Такой вывод можно сделать, проанализировав таблицу
величин погрешности 2.12.
Также стоит отметить, что в некоторых случаях есть возможность визуально по
графику изменения состояний процесса с течением времени определить наличие или от
сутствие кратных корней в знаменателе спектральной плотности стационарных процес
сов. Например, для процессов, реализация которых представлена на рисунках 2.14, 2.15,
2.16, оказалось возможно визуально определить, где присутствуют кратные корни, а где
— нет. Поэтому имеет смысл сначала проанализировать временной ряд процесса, а затем
использовать наиболее удобный метод оценки параметров и построения эффективной
модели.
44
2.5. Приложение
Рис. 2.14. Реализация процесса 2.1
Рис. 2.15. Реализация процесса 2.2
Рис. 2.16. Реализация процесса 2.3
45
Рис. 2.17. Реализация процесса 2.4
Рис. 2.18. Реализация процесса 2.5
46
Заключение
В ходе данной работы были исследованы стационарные случайные процессы с
дробно–рациональной спектральной плотностью. Также были изучены некоторые ме
тоды их моделирования.
Следует отметить, что было проведено построение моделей для процесса с дис
кретным параметром, имеющим спектральную плотность с некратными корнями зна
менателя, с корнями знаменателя второй и третьей кратности, а также с наличием
«мешающего» воздействия.
Также было проведено исследование возможности идентификации и оценки крат
ности корней знаменателя спектральной плотности стационарных моделей в дискрет
ном случае. В частности, был разработан и реализован метод, позволяющий получить
необходимые оценки для построения эффективной модели. Полученный метод был ис
следован на разнообразных примерах, что позволило определить возможности и пре
имущества разработанного алгоритма.
По проведенному исследованию можно сделать вывод о том, что в классе стаци
онарных процессов с дискретным параметром и с дробно–рациональной спектральной
плотностью есть возможность не только оценить параметры процессов и построить со
ответствующую модель, но и определить наличие кратных корней в знаменателе спек
тральной плотности и установить порядок их кратности.
Таким образом, обладая рядом наблюдений некоторого процесса, можно выбрать
наиболее удобный метод оценивания параметров и построить эффективную модель.
Это позволяет более подробно изучать и исследовать различные свойства процесса, а
также способствует проведению прогнозирования.
47
Список литературы
1. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. — М. : Наука, 1982. —
296 с.
2. Товстик Т. М. Стационарные случайные процессы с рациональными спектральными
плотностями. — СПб. : СПбГУ, 2000. — 84 с.
3. Андерсон Т. Статистический анализ временных рядов. — М. : МИР, 1976. — 760 с.
4. Бокс Д., Дженкинс Г. Анализ временных рядов. Прогноз и управление. — М. : МИР,
1974. — 406 с.
5. Носко В. П. Эконометрика. Введение в регрессионный анализ временных рядов. —
М. : Москва, 2002. — 252 с.
Отзывы:
Авторизуйтесь, чтобы оставить отзыв