close

Вход

Забыли?

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

?

637-Гайдель-3

код для вставкиСкачать
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ
Государственное образовательное учреждение высшего профессионального образования
«САМАРСКИЙ ГОСУДАРСТВЕННЫЙ АЭРОКОСМИЧЕСКИЙ УНИВЕРСИТЕТ
имени академика С. П. КОРОЛЕВА» (СГАУ)
Кафедра технической кибернетики
СТАТИСТИЧЕСКИЙ АНАЛИЗ И МОДЕЛИРОВАНИЕ ПРОЦЕССОВ АВТОРЕГРЕССИИ
И СКОЛЬЗЯЩЕГО СРЕДНЕГО
курсовая работа по дисциплине «Теория случайных процессов»
Вариант № 3
Выполнил: студент Гайдель А. В., гр. 637
Проверил: профессор Храмов А. Г.
Оценка:
Самара 2009
2
Задание
Дана реализация стационарного в широком смысле эргодического случайного
процесса с дискретным временем (стационарная случайная последовательность,
временной ряд) – выборка из 5000 последовательных значений (отсчётов) процесса.
1. Оценить моментные функции случайного процесса, рассчитав выборочное среднее,
выборочную дисперсию и выборочную нормированную корреляционную функцию.
Оценить радиус корреляции случайного процесса. Изобразить графически оценку
нормированной корреляционной функции.
2. Построить модели авторегрессии (АР), модели скользящего среднего (СС) и
смешанные модели авторегрессии и скользящего среднего (АРСС) до третьего
порядка включительно: АРСС (M, N), M = 0, 1, 2, 3; N = 0, 1, 2, 3. Каждую из
построенных моделей записать в явном виде с численными значениями параметров.
3. Рассчитать теоретические нормированные корреляционные функции выходной
последовательности для каждой из построенных выше моделей. На основе сравнения
выборочной и теоретических нормированных корреляционных функций выбрать
наиболее адекватную модель АРСС случайного процесса. Построить графики
теоретических нормированных корреляционных функций для трёх наилучших из
классов АР, СС и АРСС.
4. Построить и изобразить графически
плотности для трёх наилучших моделей.
параметрическую
оценку спектральной
5. Смоделировать случайный процесс АРСС с использованием наилучшей модели.
Сравнить графически фрагменты реализаций исходного и смоделированного
процессов.
6. Построить оценки моментных функций смоделированного процесса, сравнить их с
оценками моментных функций исходного процесса и с теоретическими моментными
функциями, соответствующими выбранной модели АРСС.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
3
Аннотация
В данной курсовой работе проводится всестороннее исследование выборки из отсчётов
некоторого неизвестного стационарного эргодического случайного процесса и
моделирование нового процесса, подобного исходному, с использованием моделей
авторегрессии и скользящего среднего различных порядков. Модели АРСС исследуются
на качество, проводится построение графиков спектральной плотности мощности для
исходного и смоделированных процессов. Для наглядности большинство результатов
изображено графически и в виде таблиц. Программа, приведённая в приложении к
курсовой работе, может служить основой для исследования любого стационарного
эргодического случайного процесса и построения моделей АРСС любых порядков.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
4
Содержание
Задание............................................................................................................................................2
Аннотация ......................................................................................................................................3
Содержание ....................................................................................................................................4
1 Оценка моментных функций ................................................................................................5
2 Построение моделей ..............................................................................................................6
3 Анализ моделей .....................................................................................................................8
4 Спектральная плотность мощности ...................................................................................13
5 Моделирование ....................................................................................................................15
6 Анализ смоделированных процессов ................................................................................15
Выводы .........................................................................................................................................17
Список использованных источников.........................................................................................18
Приложение A. Текст программы ..............................................................................................19
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
5
Оценка моментных функций
1
Пусть дана выборка x из n отсчётов стационарного в широком смысле эргодического
дискретного случайного процесса  k . Оценим моментные функции этого процесса.
Выборочное среднее – оценка математического ожидания M  k   m –
рассчитывается по формуле:
1 n
 x , где xk – соответствующие компоненты вектора x .
n k 1 k
В нашем случае, выполнив расчеты, получаем x  10.1545. Здесь и далее все данные
x
расчётов приводятся с определённой разумной степенью точности.
Формула для расчёта выборочной дисперсии имеет вид:
S2 


2
1 n
xk  x ,

n k 1
но мы воспользуемся исправленной дисперсией, рассчитывающейся по формуле



2
1 n
S2 
xk  x , так как она является лучшей (несмещённой) оценкой дисперсии

n  1 k 1
нашего процесса D k   D .

После несложных расчётов получаем: S 2  333.8049.
В соответствии с оценкой дисперсии выберем формулу для оценки корреляционной
функции R k  :

R k  



1 n k
 x  x x j k  x .
n  k  1 j 1 j
Это исправленная выборочная корреляционная функция.
В таблицу 1 запишем 11 первых неотрицательных значений этой функции, помня что


k  Z R k   R  k  :


Заметим, что R 0  S 2 , чего и следовало ожидать.
Оценим также и нормированную корреляционную функцию, которая поможет
количественно оценить корреляцию сечений:

R k 

r k   

R 0

R k 
 .
S2
Её первые значения также занесём в таблицу 1.
Таблица 1 – Первые значения выборочных корреляционных функций
0
1
2
3
4
5
k

R k  333.8049 -26.7346 -124.7274 55.6273 2.4389 1.0393

1.0000
-0.0801
-0.3737
r k 
6
7
8
k

-5.4030
9.4987
-0.7981
R k 

0.0285
-0.0024
r k  -0.0162

Заметим, r 0  1, что опять же вполне логично.
0.1666
0.0073 0.0031
9
10
-9.0512 4.0367
-0.0271 0.0121
Оценим радиус корреляции случайного процесса  k по формуле



Tcorr  min T | m  T r m  e 1 , только вместо r m используем r m , а m будем
полагать не слишком большим по сравнению с размером выборки во избежание
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
6
чрезмерной ошибки:




Tcorr  min T | m T  m  n  r m  e 1 .

В результате вычислений, полагая m  0.01 n , получаем: Tcorr  3 .
Изобразим графически на рисунке 1 оценку корреляционной функции:
Рисунок 1 – Оценка нормированной корреляционной функции
2
Построение моделей
Теперь построим все модели АРСС k 
N
M
i 0
j 1
ik i    jk  j для порядков N и M не
превышающих 3. Здесь  k – это белый шум с нулевым математическим ожиданием и
единичной дисперсией. Для простоты генерации будем считать его сечения
распределёнными по нормальному закону.
Каждый раз будем руководствоваться следующим методом:
1. Сначала отыщем коэффициенты    j j 1,M из системы линейных уравнений
 
M


k  1, M R N  k     j R N  k  j  .
j 1
2. Далее подставляем  в систему
N
M


k  0, N R k    i R i  k     j R  j  k  .
i k
j 1
Здесь R k  – смешанная корреляционная функция процессов  k и  k ,
формально выражающаяся как
t  N0 R k   M t  m  t k  m , где m и m – математические
ожидания соответствующих случайных процессов, при чём, как указано выше,
m  0 . Изначально R k  нам, вообще говоря, неизвестны.



Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
7
3. Для отыскания R k  воспользуемся следующей системой уравнений:
k
k  0, N R k    k    j R k  j  . Для определённости R 0  0 .
j 1
Подставляя R k  из этой системы в предыдущую, получаем нелинейную
систему уравнений относительно   i i0, N .
Выпишем для примера системы уравнений для моделей АР(3), СС(3) и АРСС(3, 3).
Для модели АР(3) получаем:




R 0  1 R 1   2 R 2   3 R 3   02




R 1  1 R 0   2 R 1   3 R 2
, т. е.




R 2  1 R 1   2 R 0   3 R 1




R 3  1 R 2   2 R 1   3 R 0
- 26.73461 - 124.7274 2  55.6273 3   02  333.8049

333.80491 - 26.7346 2 - 124.7274 3  -26.7346

- 26.73461  333.8049 2 - 26.7346 3  -124.7274
- 124.72741 - 26.7346 2  333.8049 3  55.6273
Для модели СС(3):

R 0   02  12   22   32

R 1  01  1 2   2 3
, т. е.

R 2  0 2  1 3

R 3  0 3
 02  12   22   32  333.8049

 01  1 2   2 3  -26.7346

 0 2  1 3  -124.7274
 0 3  55.6273
Для модели АРСС(3, 3):
R 0   0 R 0  1 R 1   2 R 2   3 R 3  1 R1   2 R2   3 R3

R 1  1 R 0   2 R 1   3 R 2  1 R0   2 R1   3 R2
R 2   R 0   R 1   R1   R0   R1
2 
3 
1
2
3
 
,
R 3   3 R 0  1 R2   2 R1   3 R0
R 4   R3   R2   R1
1
2
3
 
R 5  1 R4   2 R3   3 R2

R 6  1 R5   2 R4   3 R3
т. е.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
8
 0 R 0  1 R 1   2 R 2   3 R 3 - 26.73461 - 124.7274 2  55.6273 3  333.8049

1 R 0   2 R 1   3 R 2  333.80491 - 26.7346 2 - 124.7274 3  -26.7346
 R 0   R 1 - 26.7346  333.8049 - 26.7346  -124.7274
3 
1
2
3
 2 
,
 3 R 0 - 124.72741 - 26.7346 2  333.8049 3  55.6273
55.6273 - 124.7274 - 26.7346  2.4389
1
2
3

2.43891  55.6273 2 - 124.7274 3  1.0393

1.03931  2.4389 2  55.6273 3  -5.4030
R 0   0

R 1   0 1  1
где 
2
R 2   0 1   2  1 1   2
R 3    3  2       2       
0
1
1 2
3
1
1
2
2 1
3
 






Последняя система не была явно подставлена в основную систему по причине
нехватки места на листе, чтобы разместить получающиеся в результате такой подстановки
уравнения.
Все системы решаем численными методами. Полученные модели проверяем на
устойчивость исходя из условия, что все корни z характеристического уравнения
M

k 1
k
z k  1 лежат внутри единичной окружности z  1на комплексной плоскости.
Результаты вычислений занесём в таблицу 2.
Таблица 2 – Результат построения моделей АРСС
Параметры модели
N
M
Порядок





 jk  j
модели
k
i k i


i 0
M
0
0
0
0
1
1
1
1
2
2
2
2
3
3
3
3
3
N
0
1
2
3
0
1
2
3
0
1
2
3
0
1
2
3
0
18.2703
-1.4680
-9.7289
13.7652
18.2116
13.6406
16.8266
16.7116
-4.1868
13.6521
16.7215
14.4362
j 1
1
2
3
18.2113
-8.6484
4.4969
12.8202
-10.3812
4.0411
1
2
-0.0801
Модель существует, но не устойчива
Модель не существует
3.9211 -10.5418 4.4790
0.0438
-0.1107 -0.3825
4.9417
-0.3597 -0.4025
6.0857
16.4679
-0.4033 -0.1994
3.3721 -10.9224 4.8177
0.0781
0.0153
-0.0680 -0.3702
Модель существует, но не устойчива
Модель не существует
8.2885
-6.4464
2.5453
-0.3429 -0.1545
3
0.1116
-0.0839
Анализ моделей
Будем анализировать качество построенных моделей, сравнивая их нормированную
корреляционную функцию rMN с оценкой нормированной корреляционной функции
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
9

исходного процесса r . Для сравнения возьмём первые 10 значений нормированных
корреляционных функций и для каждой модели вычислим среднее квадратическое
отклонение
2
 MN
  rMN m  r m2 .
10

m1
2
Результаты вычислений  MN
занесём в таблицу 3.
2
Таблица 3 – Оценка качества моделей АРСС через СКО  MN
N
M
0
1
2
3
0.1758 0.1694 0.0298 0.0020
0
0.1744
0.0020
1

0.0292
0.0182
0.0040
0.0020
2
0.0262
0.0015
3

В таблице тёмно-бирюзовым цветом отмечена лучшая модель АРСС, а серым цветом –
лучшие модели из классов АР и СС. Как и следовало ожидать, это модели АРСС(3, 3),
АР(3) и СС(3) соответственно. Видно также, что все модели с СС-составляющей равной 3
оказались очень неплохими, чего нельзя сказать о соответствующих моделях с АРсоставляющей равной 3.
Теперь построим графики теоретических нормированных корреляционных функций
для указанных наилучших моделей и изобразим их на рисунках 2, 3 и 4 соответственно.
Будем считать, что для всякой модели АРСС(M, N) (N+M+1) значение нормированной
корреляционной функции совпадает, а остальные значения отыщем из системы:
M
k  N r N  M  k     j r N  M  k  j .
j 1
Рисунок 2 – Теоретическая нормированная корреляционная функция для модели АР(3)
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
10
Рисунок 3 – Теоретическая нормированная корреляционная функция для модели СС(3)
Рисунок 4 – Теоретическая нормированная корреляционная функция для модели
АРСС(3, 3)
Смоделируем три процесса на основе лучших моделей, руководствуясь формулой из
таблицы 2. Белый шум будем генерировать как последовательность из копий случайной
величины, распределённой по нормальному закону с нулевым математическим
ожиданием и единичной дисперсией. Несколько первых значений генерируемого
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
11
случайного процесса положим нулевыми как основание рекурсии. Чтобы придать новому
процессу нужное математическое ожидание (непосредственно после генерации оно
близко к нулю), просто прибавим его к каждому отсчёту. Таким образом сгенерируем
6000 отсчётов, после чего отбросим первую 1000 как брак.
Построим ещё 3 рисунка, разместив на каждом из них график выборочной
нормированной корреляционной функции исходной выборки, график параметрической
оценки нормированной корреляционной функции для соответствующей модели и график
выборочной нормированной корреляционной функции для сгенерированного по этой
модели процесса. На рисунке 5 изобразим соответствующий график для модели АР(3), на
рисунке 6 – для модели СС(3), а на рисунке 7 – для АРСС(3, 3). Чтобы 3 графика не
слишком мешали друг другу, придётся изобразить их линейную интерполяцию. При этом
следует помнить, что на самом деле нормированные корреляционные функции
дискретных процессов определены только в целых координатах, а всё остальное – это
просто их линейное продолжение, которое сделано исключительно для наглядности.
Рисунок 5 – Анализ нормированной корреляционной функции модели АР(3)
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
12
Рисунок 6 – Анализ нормированной корреляционной функции модели СС(3)
Рисунок 7 – Анализ нормированной корреляционной функции модели АРСС(3, 3)
На этих рисунках красные линии – выборочные нормированные корреляционные
функции исходной выборки, зелёные – параметрические оценки нормированных
корреляционных функций моделей, синие – выборочные нормированные корреляционные
функции соответствующих смоделированных процессов.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
13
4
Спектральная плотность мощности
Построим и изобразим на рисунках 8, 9 и 10 параметрическую оценку 
нормированной спектральной плотности мощности для трёх рассматриваемых моделей
АР(3), СС(3) и АРСС(3, 3) соответственно. Саму оценку спектральной плотности
мощности будем искать в виде:
N
  
 k eik
k 0
M
1   k e
2
.
ik
k 1
На тех же рисунках изобразим выборочную оценку нормированной спектральной
плотности мощности для исходного процесса и для соответствующего смоделированного
процесса. Такую оценку найдём из преобразования Фурье   

 Rme 
i m
, используя
m
первые 50 значений соответствующих оценок корреляционных функций вместо
бесконечного числа. Все оценки спектральной плотности мощности нормируются на
соответствующие оценки дисперсий.
Рисунок 8 – Спектральная плотность мощности для модели АР(3)
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
14
Рисунок 9 – Спектральная плотность мощности для модели СС(3)
Рисунок 10 – Спектральная плотность мощности для модели АРСС(3, 3)
Цветографическая схема та же, что и на рисунках 5 – 7: красным отмечена оценка
нормированной спектральной плотности мощности для исходного процесса, зелёным –
параметрическая оценка для модели, а синим – выборочная оценка для соответствующего
смоделированного процесса.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
15
5
Моделирование
Итак, лучшая модель получилась АРСС(3, 3). Приведём на рисунке 8 фрагменты
реализаций исходного и построенного процесса из 100 первых отсчётов (не бракованных,
конечно).
Рисунок 11 – Фрагмент реализации случайного процесса модели АРСС(3, 3)
На этом рисунке красным цветом показана реализация исходного процесса, синим –
смоделированного. Кроме того, для наглядности голубым цветом отмечено выборочное

среднее x исходного процесса, а зелёным – прямые вида x  S . Видно, что характер
реализаций совпадает.

6

Анализ смоделированных процессов
Наконец, приведём итоговую таблицу 4, в которую соберём все статистические и
теоретические сведения, необходимые для наглядного анализа трёх лучших моделей.
Таблица 4 – Итоговый анализ построенных моделей
Параметры Исходный
АР(3)
СС(3)
АРСС(3, 3)
процесса
процесс
Теория Выборка Теория Выборка Теория Выборка
-45.3730
-65.3317
-56.5194
-64.7280
Минимум
67.1167
72.4522
75.2665
Максимум 68.3770
10.1545
10.1545 10.1294
10.1545 10.0735
10.1545 9.9608
Среднее
333.804 321.0027 333.8049 336.8887 333.8049 336.8630
Дисперсия 333.8049
9
18.2703 17.9165
18.2703 18.3545
18.2703 18.3538
Стандартное 18.2703
отклонение
Нормированная корреляционная функция
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
r(0)
-0.0801
-0.0801 -0.0606
-0.0801
-0.0855
-0.0801
-0.0956
r(1)
-0.3737
-0.3737 -0.3717
-0.3737
-0.3908
-0.3737
-0.3873
r(2)
0.1666
0.1666
0.1479
0.1666
0.1877
0.1666
0.1799
r(3)
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
16
r(4)
r(5)
r(6)
r(7)
r(8)
r(9)
r(10)
СКО
0.0073
0.0031
-0.0162
0.0285
-0.0024
-0.0271
0.0121
0.0000
0.1180
-0.1114
-0.0175
0.0556
-0.0097
-0.0219
0.0113
0.0262
0.1110
-0.1048
-0.0118
0.0489
-0.0208
-0.0167
0.0226
0.0241
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0020
0.0130
-0.0226
0.0097
0.0090
0.0069
-0.0070
-0.0142
0.0037
0.0073
0.0031
-0.0162
0.0045
0.0007
0.0004
-0.0006
0.0015
0.0039
-0.0044
-0.0061
0.0028
-0.0050
0.0065
0.0057
0.0026
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
17
Выводы
Задачи моделирования случайных процессов возникают на практике довольно часто.
Это, прежде всего, связано с экономикой её экономическими процессами. Модели
авторегрессии и скользящего среднего позволяют моделировать случайные процессы,
подобные исходному, по уже имеющейся реализации такого исходного процесса. Общая
модель, предложенная Боксом и Дженкинсом, включает как параметры авторегрессии, так
и параметры скользящего среднего. Она позволяет добиться максимального подобия
новых смоделированных процессов.
В этой работе было проведено исследование подобного рода моделирования для
некоторого исходного неизвестного эргодического процесса. В ходе работы была
проанализирована выборка из отсчётов исходного процесса, построены все смешанные
модели АРСС до третьего порядка включительно, проведён поиск наилучшей модели и
моделирование нового случайного процесса по ней. При этом написана универсальная
программа, позволяющая строить, вообще говоря, смешанные модели любых порядков.
Методы исследования, использующиеся в работе, могут быть применены на практике
для реального статистического анализа и моделирования любого эргодического
случайного процесса. Подобные методы моделирования актуальны на сегодняшний день и
находятся на стадии исследования.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
18
Список использованных источников
1. Тараскин, А.Ф. Статистический анализ временных рядов авторегрессии и
скользящего среднего: учебное пособие [Текст] // Самара: СГАУ, 1998. – 56с.
2. Тараскин, А.Ф. Статистическое моделирование и метод Монте–Карло: учебное
пособие [Текст] // Самара: СГАУ, 1997. – 62с.
3. Храмов, А.Г. Анализ и моделирование процессов АРСС: интернет-ресурс к курсовой
работе [Электронный ресурс] // Самара: СГАУ, 2009.
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
19
Приложение A. Текст программы
clear();
// Initial settings
SOURCE_FILE_NAME = 'D:\temp\source.txt'; // File contains input data
FLOAT_FORMAT = '%16.4f'; // Floating point values representation format
INT_FORMAT = '%d'; // Decimal integer values representation format
EPSILON = 1.0E-6; // Precision
MAX_AR_LEVEL = 3;
MAX_MA_LEVEL = 3;
IMITATION_LENGTH = 5000;
// Input data
x = fscanfMat(SOURCE_FILE_NAME); // Sample
n = length(x); // Sample size
// Helpers
function printMat(M, mformat),
[n, m] = size(M);
for i = 1:n,
for j = 1:m,
printf(mformat + " ", M(i, j));
end;
printf("\n");
end;
endfunction;
// 1. Moment functions
// Correlation function esimation
function R = correlation(k, x)
if (k < 0),
k = -k;
end;
R = 0.0;
n = length(x);
meanx = mean(x);
for i = 1:(n-k),
R = R + (x(i) - meanx) * (x(i+k) - meanx);
end;
R = R / (n - k - 1);
endfunction;
// Normalized correlation function
function r = ncorrelation(k, x)
r = correlation(k, x) / correlation(0, x);
endfunction;
// Correlation distance
function T = corrdist(x)
coefficient = 0.01;
T = coefficient * length(x) - 2;
em1 = exp(-1);
while (T >= 0) & (abs(ncorrelation(T, x)) < em1),
T = T - 1;
end;
T = T + 1;
endfunction;
// Drawing normalized correlation function plot
function corrplot()
p.thickness = 6;
m = 10;
t = [-m:m];
y = zeros(length(t), 1);
for i = 1:length(t),
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
20
y(i) = ncorrelation(t(i), x);
end;
plot2d3(t, y, axesflag=5, style=2);
a = gca();
p = a.children.children;
p.thickness = 3;
p.mark_mode = "on";
p.mark_size_unit = "point";
p.mark_style = 11;
p.mark_size = 3;
endfunction;
// Main
meanx = mean(x); // Sample mean
svx = variance(x); // Sample variance
m = 10; // Correlation values count
R = zeros(m, 1);
r = zeros(m, 1);
for i = 0:m,
R(i+1) = correlation(i, x);
r(i+1) = ncorrelation(i, x);
end;
Tcorr = corrdist(x);
printf("Sample mean: " + FLOAT_FORMAT + "\n", meanx);
printf("Sample variance: " + FLOAT_FORMAT + "\n", svx);
printf("Correlation function estimation:\n");
printMat(R, FLOAT_FORMAT);
printf("Normalized correlation function estimation:\n");
printMat(r, FLOAT_FORMAT);
printf("Correlation distance: " + INT_FORMAT + "\n", Tcorr);
scf(1);
corrplot();
// 2. Models building
// Autoregression coefficients search
function betas = ar(x, arLevel, maLevel)
R = zeros(2 * arLevel, 1);
for i = (maLevel - arLevel + 1) : (maLevel + arLevel),
R(i - maLevel + arLevel) = correlation(i, x);
end;
Rmm = zeros(arLevel, arLevel);
_n = maLevel;
for i = 1:arLevel,
_m = _n;
for j = 1:arLevel,
Rmm(i, j) = R(_m - maLevel + arLevel);
_m = _m - 1;
end;
_n = _n + 1;
end;
Rm = -R(arLevel + 1 : 2 * arLevel);
betas = linsolve(Rmm, Rm);
endfunction;
// Mutual Correlation Function
function rrr = mcorrelation(k, alph, betas)
rrr = alph(k+1);
len = min(k, length(betas));
for j = 1 : len,
rrr = rrr + betas(j) * mcorrelation(k - j, alph);
end;
endfunction;
// Moving average coefficients search
function alphas = ma(x, arLevel, maLevel, betas)
for i = 0 : max([arLevel, maLevel]),
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
21
R(i+1) = correlation(i, x);
end;
function zr = syst(alph)
for k = 0 : maLevel,
zr(k+1) = -R(k+1);
for i = k : maLevel,
zr(k+1) = zr(k+1) + alph(i+1) * mcorrelation(i - k, alph, betas);
end;
for j = 1 : arLevel,
zr(k+1) = zr(k+1) + betas(j) * R(abs(k - j) + 1);
end;
end;
endfunction;
[alphas, values, info] = fsolve([1 : (maLevel+1)], syst);
for i = 1:length(values),
if (abs(values(i)) > EPSILON | info == 4) then
alphas(1) = %i;
break;
end;
end;
endfunction;
// Image vector
function s = image(v)
s = %F;
for i = 1:length(v),
if (imag(v(i)) <> 0) then
s = %T;
break;
end;
end;
endfunction;
// Model stability
function s = stable(betas)
p = poly([pertrans(-betas) 1], "z", "coeff");
z = roots(p);
s = %T;
for i = 1:length(z),
if (abs(z(i)) >= 1) then
s = %F;
break;
end;
end;
endfunction;
// Main
alphas_list = list();
betas_list = list();
for i = 0 : MAX_AR_LEVEL,
for j = 0 : MAX_MA_LEVEL,
betas = ar(x, i, j);
alphas = ma(x, i, j, betas);
alphas_list($+1) = alphas;
betas_list($+1) = betas;
printf("ARMA(" + INT_FORMAT + "," + INT_FORMAT + ")\n", i, j);
if (image(alphas)) then
printf("Model does not exist.\n");
continue;
end;
if (~stable(betas)) then
printf("Model exists, but not stable.\n");
continue;
end;
printf("alpha:\n")
printMat(alphas, FLOAT_FORMAT);
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
22
printf("beta:\n");
printMat(betas, FLOAT_FORMAT);
end;
end;
// 3. Models Analysis
// Theoretical correlation function for ARMA model
function R = theoretical_corr(betas, startR, k)
nm = length(startR) - 1;
k = abs(k);
if (k > nm) then
R = 0;
M = length(betas);
for j = 1 : M,
R = R + betas(j) * theoretical_corr(betas, startR, k - j);
end;
else
R = startR(k + 1);
end;
endfunction;
// Normalized theoretical correlation function for ARMA model
function r = norm_theoretical_corr(betas, startR, k)
r = theoretical_corr(betas, startR, k) / theoretical_corr(betas, startR,
0);
endfunction;
// Quadratic Error
function epsilon = quadratic_error(x, y)
epsilon = 0;
m = min(length(x), length(y));
for j = 1 : m,
epsilon = epsilon + (x(j) - y(j))^2;
end;
endfunction;
function corrplot2(betas, R, N, M, m_)
p.thickness = 6;
m = m_;
t = [-m:m];
y = zeros(length(t), 1);
for i = 1:length(t),
y(i) = norm_theoretical_corr(betas, R(1 : N + M + 1), t(i));
end;
plot2d3(t, y, axesflag=5, style=2);
a = gca();
p = a.children.children;
p.thickness = 3;
p.mark_mode = "on";
p.mark_size_unit = "point";
p.mark_style = 11;
p.mark_size = 3;
endfunction;
function eta = imitate(alphas, betas, meanx, count)
defect = 1000;
eta = zeros(count + defect + 1, 1);
ksi = grand(count + defect + 1, 1, 'nor', 0, 1);
N = length(alphas) - 1;
M = length(betas);
for k = 1 : count + defect + 1,
eta(k) = 0;
for i = 0 : N,
if (k - i > 0) then
eta(k) = eta(k) + alphas(i+1) * ksi(k - i);
end;
end;
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
23
for j = 1 : M,
if (k - j > 0) then
eta(k) = eta(k) + betas(j) * eta(k - j);
end;
end;
end;
eta = eta(defect + 2 : count + defect + 1) + meanx;
endfunction;
m = 10; // Analysis Depth
epsilon = zeros(MAX_AR_LEVEL, MAX_MA_LEVEL);
best_ar_eps = %inf;
best_ma_eps = %inf;
best_arma_eps = %inf;
best_ar_alpha = [];
best_ma_alpha = [];
best_arma_alpha = [];
best_ar_beta = [];
best_arma_beta = [];
R = zeros(m+1, 1);
r = zeros(m+1, 1);
for k = 0 : m,
R(k+1) = correlation(k, x);
r(k+1) = R(k+1) / R(1);
end;
r_model = zeros(m+1, 1);
for i = 0 : MAX_AR_LEVEL,
for j = 0 : MAX_MA_LEVEL,
alphas = alphas_list(i * (MAX_MA_LEVEL + 1) + j + 1);
betas = betas_list(i * (MAX_MA_LEVEL + 1) + j + 1);
for k = 0 : m,
r_model(k+1) = norm_theoretical_corr(betas, R(1 : (i + j + 1)), k);
end;
epsilon(i+1, j+1) = quadratic_error(r, r_model);
if (i == 0) & (epsilon(1, j+1) < best_ma_eps) then
best_ma_eps = epsilon(1, j+1);
best_ma_alpha = alphas;
elseif (j == 0) & (epsilon(i+1, 1) < best_ar_eps) then
best_ar_eps = epsilon(i+1, 1);
best_ar_alpha = alphas;
best_ar_beta = betas;
elseif (epsilon(i+1, j+1) < best_arma_eps) then
best_arma_eps = epsilon(i+1, j+1);
best_arma_alpha = alphas;
best_arma_beta = betas;
end;
end;
end;
printf("Epsilon:\n");
printMat(epsilon, FLOAT_FORMAT);
printf("Best models:\nAR(" + INT_FORMAT + "), MA(" + INT_FORMAT + "), ARMA("
+ INT_FORMAT + "," + INT_FORMAT + ").\n", length(best_ar_beta),
length(best_ma_alpha) - 1, length(best_arma_beta), length(best_arma_alpha) 1);
scf(2);
corrplot2(best_ar_beta, R, length(best_ar_beta), 0, m);
scf(3);
corrplot2([], R, 0, length(best_ma_alpha) - 1, m);
scf(4);
corrplot2(best_arma_beta, R, length(best_arma_beta), length(best_arma_alpha)
- 1, m);
eta_ar = imitate(best_ar_alpha, best_ar_beta, meanx, IMITATION_LENGTH);
eta_ma = imitate(best_ma_alpha, [], meanx, IMITATION_LENGTH);
eta_arma = imitate(best_arma_alpha, best_arma_beta, meanx, IMITATION_LENGTH);
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
24
r_ar = zeros(m+1, 1);
r_ma = zeros(m+1, 1);
r_arma = zeros(m+1, 1);
r_ar_imit = zeros(m+1, 1);
r_ma_imit = zeros(m+1, 1);
r_arma_imit = zeros(m+1, 1);
for k = 0 : m,
r_ar(k + 1) = norm_theoretical_corr(best_ar_beta, R(1 :
length(best_ar_beta) + 1), k);
r_ma(k + 1) = norm_theoretical_corr([], R(1 : length(best_ma_alpha)), k);
r_arma(k + 1) = norm_theoretical_corr(best_arma_beta, R(1 :
length(best_arma_alpha) + length(best_arma_beta)), k);
r_ar_imit(k + 1) = ncorrelation(k, eta_ar);
r_ma_imit(k + 1) = ncorrelation(k, eta_ma);
r_arma_imit(k + 1) = ncorrelation(k, eta_arma);
end;
scf(5);
plot2d([0 : m], [r r_ar r_ar_imit], style=[5 3 2], axesflag=5,
leg="Source@AR(" + string(length(best_ar_beta)) + ")@Imitation");
scf(6);
plot2d([0 : m], [r r_ma r_ma_imit], style=[5 3 2], axesflag=5,
leg="Source@MA(" + string(length(best_ma_alpha) - 1) + ")@Imitation");
scf(7);
plot2d([0 : m], [r r_arma r_arma_imit], style=[5 3 2], axesflag=5,
leg="Source@ARMA(" + string(length(best_arma_beta)) + "," +
string(length(best_arma_alpha) - 1) + ")@Imitation");
// 4. Power Spectral Density
// Power Spectral Density
function Fi = pow_spec_dens_arma(omega, alphas, betas)
s_up = 0;
for k = 0 : length(alphas) - 1,
s_up = s_up + alphas(k+1) * exp(%i * k * omega);
end;
s_down = 1;
for k = 1 : length(betas),
s_down = s_down - betas(k) * exp(%i * k * omega);
end;
Fi = abs(s_up / s_down)^2;
endfunction;
function Fi = pow_spec_dens(omega, R)
Fi = R(1);
for k = 1 : length(R) - 1,
Fi = Fi + 2 * R(k + 1) * cos(omega * k);
end;
endfunction;
// Power Spectral Density Plot
function densplot(alphas, betas, R, R_imit)
omega = [0 : 0.01 : %pi];
len = length(omega);
dens = zeros(len, 1);
dens_source = zeros(len, 1);
dens_imit = zeros(len, 1);
for j = 1 : len,
dens(j) = pow_spec_dens_arma(omega(j), alphas, betas) / R(1);
dens_source(j) = pow_spec_dens(omega(j), R) / R(1);
dens_imit(j) = pow_spec_dens(omega(j), R_imit) / R_imit(1);
end;
str = "";
if (length(alphas) == 1) then
str = "AR(" + string(length(betas)) + ")";
elseif (length(betas) == 0) then
str = "MA(" + string(length(alphas) - 1) + ")";
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
25
else
str = "ARMA(" + string(length(betas)) + "," + string(length(alphas) - 1)
+ ")";
end;
plot2d(omega, [dens_source dens dens_imit], style=[5 3 2], axesflag=5,
leg="Source@" + str + "@Imitation");
endfunction;
perc = 0.01;
R_50 = zeros(perc * length(x), 1);
R_ar_imit_50 = zeros(perc * length(x), 1);
R_ma_imit_50 = zeros(perc * length(x), 1);
R_arma_imit_50 = zeros(perc * length(x), 1);
for k = 0 : length(R_50) - 1,
R_50(k+1) = correlation(k, x);
R_ar_imit_50(k+1) = correlation(k, eta_ar);
R_ma_imit_50(k+1) = correlation(k, eta_ma);
R_arma_imit_50(k+1) = correlation(k, eta_arma);
end;
scf(8);
densplot(best_ar_alpha, best_ar_beta, R_50, R_ar_imit_50);
scf(9);
densplot(best_ma_alpha, [], R_50, R_ma_imit_50);
scf(10);
densplot(best_arma_alpha, best_arma_beta, R_50, R_arma_imit_50);
// 5. Imitation
function imitation_plot(x, imitation, meanx, sv)
t = [1 : 100];
q = [x(t) imitation(t) (zeros(length(t), 1) + meanx) (zeros(length(t), 1) +
meanx + sqrt(sv)) (zeros(length(t), 1) + meanx - sqrt(sv))];
plot2d(t, q, style=[5 2 4 3 3], axesflag=5,
leg="Source@Imitation@Mean@Standard deviation");
endfunction;
best_alphas = [];
best_betas = [];
if (best_ar_eps < best_ma_eps) then
if (best_ar_eps < best_arma_eps) then
best_alphas = best_ar_alpha;
best_betas = best_ar_beta;
else
best_alphas = best_arma_alpha;
best_betas = best_arma_beta;
end;
else
if (best_ma_eps < best_arma_eps) then
best_alphas = best_ma_alpha;
best_betas = best_ma_beta;
else
best_alphas = best_arma_alpha;
best_betas = best_arma_beta;
end;
end;
printf("Best model: ARMA(" + INT_FORMAT + "," + INT_FORMAT + ")\n",
length(best_betas), length(best_alphas) - 1);
for i = 1 : length(r_model),
r_model(i) = norm_theoretical_corr(best_betas, R(1 : (length(best_betas) +
length(best_alphas))), i-1);
end;
imitation = imitate(best_alphas, best_betas, meanx, IMITATION_LENGTH);
scf(11);
imitation_plot(x, imitation, meanx, svx);
// 6. Imitation Analysis
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
26
function total_sample_analysis(x)
m = 10;
R = zeros(m + 1, 1);
R_src = zeros(m + 1, 1);
for k = 0 : m,
R(k + 1) = correlation(k, x);
end;
printf("Minimum: " + FLOAT_FORMAT + "\n", min(x));
printf("Maximum: " + FLOAT_FORMAT + "\n", max(x));
printf("Mean: " + FLOAT_FORMAT + "\n", mean(x));
printf("Variance: " + FLOAT_FORMAT + "\n", R(1));
printf("Standard devation: " + FLOAT_FORMAT + "\n", sqrt(R(1)));
printf("Normalized correlation function:\n");
printMat(R / R(1), FLOAT_FORMAT);
printf("Epsilon: " + FLOAT_FORMAT + "\n", quadratic_error(R / R(1), r));
endfunction;
function total_model_analysis(alphas, betas, R)
R_model = zeros(length(R), 1);
for k = 0 : m,
R_model(k + 1) = theoretical_corr(betas, R(1 : length(alphas) +
length(betas)), k);
end;
printf("Mean: " + FLOAT_FORMAT + "\n", meanx);
printf("Variance: " + FLOAT_FORMAT + "\n", R_model(1));
printf("Standard devation: " + FLOAT_FORMAT + "\n", sqrt(R_model(1)));
printf("Normalized correlation function:\n");
printMat(R_model / R_model(1), FLOAT_FORMAT);
printf("Epsilon: " + FLOAT_FORMAT + "\n", quadratic_error(R_model /
R_model(1), r));
endfunction;
mean_imit = mean(imitation);
var_imit = variance(imitation);
m = 10;
r_imit = zeros(m, 1);
for k = 0 : m,
r_imit(k+1) = ncorrelation(k, imitation);
end;
printf("Imitation mean: " + FLOAT_FORMAT + "\n", mean_imit);
printf("Imitation variance: " + FLOAT_FORMAT + "\n", var_imit);
printf("Imitation normalized correlation function:\n");
printMat(r_imit, FLOAT_FORMAT);
printf("ARMA model theoretical correlation function:\n");
printMat(r_model, FLOAT_FORMAT);
printf("\nTotal table\n");
printf("\nSource\n");
total_sample_analysis(x);
printf("\nAR(" + INT_FORMAT + ")\n", length(best_ar_beta));
printf("Theory\n");
total_model_analysis(best_ar_alpha, best_ar_beta, R);
printf("Sample\n");
total_sample_analysis(eta_ar);
printf("\nMA(" + INT_FORMAT + ")\n", length(best_ma_alpha) - 1);
printf("Theory\n");
total_model_analysis(best_ma_alpha, [], R);
printf("Sample\n");
total_sample_analysis(eta_ma);
printf("\nARMA(" + INT_FORMAT + "," + INT_FORMAT + ")\n",
length(best_arma_beta), length(best_arma_alpha) - 1);
printf("Theory\n");
total_model_analysis(best_arma_alpha, best_arma_beta, R);
printf("Sample\n");
total_sample_analysis(eta_arma);
Статистический анализ и моделирование процессов авторегрессии и скользящего
среднего. Гайдель А. В., гр. 637
Документ
Категория
Экономико-математическое моделирование
Просмотров
88
Размер файла
519 Кб
Теги
гайдель, 637
1/--страниц
Пожаловаться на содержимое документа