© Семёнычев В. К., Сергеев А. В., 2010
СЕМЁНЫЧЕВ Валерий Константинович
Доктор экономических наук, доктор технических наук, профессор, ректор
Самарский муниципальный институт управления
443084, РФ, г. Самара, ул. Стара Загора, 96 Контактный телефон: (846) 951-54-66 e-mail: 505tot@smim.ru
СЕРГЕЕВ Алексей Викторович
Аспирант кафедры математических методов и информационных технологий
Самарский муниципальный институт управления
443084, РФ, г. Самара, ул. Стара Загора, 96 Контактный телефон: (906) 125-82-65 e-mail: alexey.sergeev.samara@gmail.com
Метод МНК-оценок параметров тренд-сезонного ряда
Ключевые слова: ДКМД-модели; идентификация; статистические свойства оценок; моделирование; прогнозирование; гармоника.
Аннотация. В статье рассматриваются статистические свойства оценок параметров моделей суммы линейного тренда, колебательной компоненты и ошибки, полученных на основе параметрических моделей авторегрессии - скользящего среднего. В качестве колебательной компоненты рассматриваются одна гармоника и сумма двух гармоник. В статье получены аналитические выражения смещений оценок частот и их дисперсий. Предложен метод компенсации погрешности от смещения оценок и уменьшения их дисперсии. Проведены численные исследования предложенного метода на тестовых выборках и приведен результат применения методов на реальной выборке.
Структура рядов экономических показателей зачастую предполагается в виде суммы тренда, колебательной и стохастической компонент [1]. При реализации моделирования и прогнозирования таких тренд-сезонных рядов на относительно коротких выборках нелинейный тренд может быть представлен асимптотически в виде линейной функции (т. е. первыми двумя членами разложения в ряд Тейлора нелинейного в общем случае тренда в малой окрестности точки разложения), а колебательная компонента -гармоникой или суммой гармоник [2].
1-й случай. Идентификация модели в виде суммы линейного тренда, гармоники и ошибки
На практике широкое применение получила модель
ук = а0 + а1к + а28Іп(юк +ф) + єк, (1)
где ук - уровни временного ряда; а0, а1, а2, ш, ф - параметры модели; к - номер наблюдения (к = 1, N, N - объем выборки); £к - стохастическая компонента, для которой обычно принимают справедливыми условия Гаусса-Маркова: а) М[єк] = 0 для к = 1, N; б) для і Ф. М[є.є.] = 0, М[є2] = ст2 для к = 1, N, где М - оператор математического ожидания.
Параметрическую идентификацию модели (1) будем производить, в соответствии с [3], в два этапа: на первом - найдем оценку частоты ш гармоники; на втором - оценку
остальных параметров модели: а0, а , а , ф. Для (1) (при к > 4) справедлива параметрическая модель авторегрессии - скользящего среднего (АВМА-модель) четвертого порядка [3]:
yk - (X1 + 2)yk-1 + 2(X1 + 1yk-2 - (X1 + 2)Уk-3 + yk-4 = = Єk - (X1 + 2)ЄИ + 2(X1 + ^k-2 - (X1 + 2^k-3 +Єk-4,
(2)
где а. = cos о.
В качестве оператора сглаживания для идентификации параметров модели (2) целесообразно использовать метод наименьших квадратов (МНК) как наиболее простой и распространенный, позволяющий, при выполнении условий Гаусса-Маркова для стохастической компоненты, получить оптимальные (несмещенные, эффективные и состоятельные) оценки параметров. Обозначим
5k = Єk - (X1 + 2)Єk-1 + 2(X1 + 1)Єk-2 - (X1 + 2)Єk-3 +%-4
и применим МНК к модели:
yk - (X1 + 2)yk-1 + 2(X1 + 1yk-2 - (X1 + 2)yk-3 + yk-4 =5k .
(з)
(4)
Определим свойства МНК-оценки X* для модели (4). Из предположения а) для ек очевидно, что стохастическая компонента имеет, так же как и ек, нулевое математическое ожидание:
М[^к ] = М[ек ]- (Х1 + 2)МК-1 ]+2(Х1 + 1)м[ек-2 ]- (Х1 + 2)м[бк-з ] +
+М[&к-з ] = 0 — (X + 2)0 + 2(Х +1)0 — (X + 2)0 + 0 = 0. (5)
Ковариационная матрица для ^к имеет вид:
cov(5k) =
где “1 = 14а2 + 5А1а2 + X; “4 = -4а2 -2А1а2, “5 = а2.
11 ( “2 “3 “4 “5 0 ... 0I
0
(2 1 2 “3 4 5 ...
0
“3 2 1 2 “3 4 ...
(4 0
“3 2 1 2 “3 ...
( 0 0
|_|5 (( 0 4 “5 “3 “4 2 “3 1 “2 2 ... “1 ...
0 V 0 0 0 0 0 ... “ У
2 -12а2 4 -1 1а2 - 4X2 а2, “3
(б)
= 8а + 8^а +Х1а ,
При i < j - 4, М[є£.] Ф 0, т. е. существует автоковариация £k. Гомоскедастичность
ошибки £к (предположение Ь) обеспечивает гомоскедастичность ^к:
МКП = М[(^ — (Х1 + 2)^—1 + 2(Х1 + 1)е^ 2 — (Х1 + 2)е(—3 +е(—4)2] =
= М |^( е2 + 4е2— 1 + 4е2—2 + 4е2—3 + е2—4 +Х162—1 + 4Х1е2—2 + А1е2— 3 )2 ] = 14ст2 + 5Х1ст2 + Х2 а2. (7)
Оператор МНК для модели (4) примет вид:
X* = argmin ( —Ц-^ - (Xl + 2)yk-1 + 2(Xl +1)yk-2 - (Xl + 2)yk-з + yk-4)2 |. а. IN - 4 k=5 у
(8)
Введем обозначения:
ак = Ук - 2( У к-1 - У к-2 + Ук-з) + Ук-4> Ьк = Ук-1 - 2 У к-2 + Ук-3>
Ск = 8к - 2(ек-1 - ек-2 + 8к-3 ) + ек-4> йк =8к-1 - 2ек-2 +8к-3'
При (9-10) выражение (8) примет вид:
X* = а^тт| —Ц-Е («к -^1Ьк )2 I
V N - 4 к=5 у
(9)
(10)
(11)
(12)
(13)
Осуществляя дифференцирование по Х1
й Г 1
йXl V N - 4
получим нормальное уравнение
2
2
Л I *т лЕ(ак Х1Ьк) 1“ , . |Е акЬк Х1ЕЬ
N - 4
N - 4
откуда МНК-оценка X* будет равна
ЕакЬк Х1 ЕЬк | = 0>
к=5 к=5
X* =-
ЕаА
Еь
(14)
Для исследования смещенности оценки X* подставим в (14) выражение ак в явном виде из исходной модели:
Е(Ск Х1йк +Х1Ьк )Ьк
X* = ^5
(15)
Еь
Тогда математическое ожидание (15) примет вид: М^*] = М
ЕС (ск X1йk +X1bk )Ьк Е (ск ^'"Ак )Ьк 1 1
к=5 =М к=5 + X1М к=5
I 2к ♦■О 1 ЕЬ2 _ к=5 _ ЕЬ2 _ к=5 _
= М
Е(ск -Мк )Ьк
Еь
+ Xl'
к=5
к=5
к=5
к=5
к=5
к=5
к=5
к=5
Видим, что оценка X* смещена на величину
N
Х(ск -ХА Ж
м
Хж
(16)
При этом величина Ьк в (10) будет мала при относительно малом отличии друг от друга соседних четырех уровней ряда и, наоборот, будет велика при значительном их отличии. В отношении стохастической компоненты то же можно сказать о величинах ск в (11) и йк в (12) при соседних четырех и пяти уровнях ряда соответственно. Смещение (16) можно уменьшать путем увеличения Ьк и уменьшения ск и йк, т. е. необходимо достичь значительного отличия уровней ряда и малого отличия стохастической компоненты. Для этого целесообразно применить прием прорежения выборки: удаления из рассмотрения (расчета) каждого г-го наблюдения, в результате чего получатся г прореженных выборок. Одновременно необходимо накладывать на полученную таким образом выборку условие равноудаленности соседних наблюдений. На каждой полученной выборке производится идентификация Х*г выражением (14). Если прорежение производится таким образом, что для расчета оценки частоты остается каждое А наблюдение (шаг прорежения), то оценку рассчитываем по формуле
Х1 = 2со8(юД).
(17)
Из идентифицированных на каждой из прореженных таким образом г-й выборке оценок Х*г целесообразно выбирать ту, которая обеспечивает минимум остаточной дисперсии. Прорежение с шагом А = 5 обращает автоковариационную матрицу в матрицу с нулевыми элементами и ненулевыми диагональными элементами, равными дисперсии ^к. Нетрудно показать, что прорежение выборки уменьшает дисперсию оценки X*; (В[] - оператор дисперсии):
В[Х*] = В
Х(Ск -ХА + Х1Жк)Жк
=м
( N
Х(Ск -ХА + Х1Жк)Жк
^--------*-------------м
ХК
5
N
Х(ск -ХА + Х1Жк)Жк
Хж
ХЖ
Х (ск -ХА )Жк Х (ск -ХА )Жк Х (ск -ХА )Ьк
+Х1 - м к=5 N -Х1 = В к=5 N
ХК _ к=5 ХК _ к=5 _ ХК _ к=5 _
=м
Оставшиеся параметры модели (1) а0, а , а2, ф оцениваются МНК при известной оценке X* и стохастической компоненте ек из системы линейных алгебраических уравнений (СЛАУ), определяемой условием
й0> а1> й3’ а4 ^ к = 1
а0, а*, а*, а* = аг§тт<! Х(Ук - а0 - а1к - а* зт(ю*к) - я4 соз(ю*к)) к
(18)
к=5
к=5
к=5
2
к=5
к=5
к=5
где аз = Я2С08 (ф), Я4 = Я^Ш (ф). Отсюда:
= ^(я* )2 + (я* )2,
ф = ягс^
'О
+ ли, и е X.
(19)
Статистические свойства оценок я0, я1, я2, ф в (18), (19) известны и широко рассмотрены в современной эконометрической литературе, что нельзя сказать об оценках, получаемых на первом этапе идентификации.
В качестве примера приведем результаты моделирования уровня безработицы моделью суммы линейного тренда и одной гармоники с использованием прорежения и без него (рис. 1).
Номер наблюдения
--- Количество не занятых трудовой деятельностью
в Самарской области (2004-2006 гг.)
... Моделирование без прорежения (коэффициент детерминации 0,73)
---Моделирование с прорежением (коэффициент детерминации 0,96)
Рис. 1. Моделирование моделью суммы линейного тренда и одной гармоники с прорежением и без прорежения выборки
2-й случай. Идентификация модели в виде суммы линейного тренда, двух гармоник и ошибки
Для модели
ук = я0 + я1к + Х я1^ш^к + ф() + ек
1=1
справедлива при к > 7 параметрическая АВ.МА-модель шестого порядка:
ук - (2 + Х1 +Х2)ук-1 + (31 + 2(Х1 +Х2) +Х1Х2)Ук—2 --(4 + 2(^1 +^2) + Х1^2)Ук-з + (3 + 2(^1 +^2) + ^1^2)Ук—4 -
-(2 + Я,1 + ^2)ук-5 + ук-6 = ек - (2 + Х1 + Х2)8к-1 +
+(31 + 2(Х1 +Х2) +Х1Х2)8к-2 - (4 + 2(Х1 +Х2) +Х1Х2)8к-3 + +(3 + 2(^1 +^2) +Х1Х2)8к-4 - (2 + Х1 +Х2)8к-5 +8k-6,
(20)
где Х1 = 2с08 (ш1), Х2 = 2с08 (ш2).
V яз У
Обозначим скользящее среднее в (21) через стохастическую компоненту ^к:
5к =8к - (2 + Х1 +Х2)8к-1 + (31 + 2(Х1 +Х2) + Х1Х2)8к-2 -
-(4 + 2(Х1 +Х2) + Х1Х2)ек-3 +(3 + 2(Х1 + Х2) + Х1Х2)ек-4 -(2 + Х1 +Х2)ек-5 + ек-6. (22)
Относительно ошибки ^к можно сделать выводы о центрированности, существовании ненулевой автоковариации ^к при соседних шести наблюдениях и гомоскедастич-ности. Введем следующие обозначения:
т = Х1 + Х2, (23)
g = ^ (24)
Як = ук - 2(ук-1 + ук-5) + 3(ук-2 + ук-4) - 4(ук-3 + Я-Л
Ьк = -ук-1 + 2(ук-2 - ук-3 + ук-4) - ук—5 ’
Ск = у к-2 - 2ук-3 + ук-4’
Ак = £к - 2(£к-1 + £к-5) + 3(£к-2 + £к-4) - 4£к-3 + £к-6, ек = -£к-1 + 2(£к-2 - £к-3 + £к-4) - £к-5’ fk = £к-2 - £к-3 + £к-4.
Оператор МНК для модели (22) во введенных обозначениях примет вид:
т\ g * = а^шт | Т7~7У(як + тЬк + gck )2 |. (25)
т, g V ^ - 6 к=7 )
Полученные оценки т и g подставляем в уравнение х2 - тх - g = 0. Если корни вещественные и принадлежат -2 < т < 2, -2 < g < 2, то решение будет найдено из системы уравнений:
Гт = Х1 + Х2,
1 g = Х1Х2.
Иначе, решение лежит на границе квадрата. То есть необходимо рассмотреть 2 случая, так как Х1 и Х2 входят в выражение симметрично.
N 2
.2 = а^шт У (як + (2 + ^2 )Ьк + 2^2Ск'
Если X, = 2, то х2 = агБш1п У (як + (2 + Х2 Ж + 2Х2Ск)
Х2 к=7
N
Если Х1 = -2, то = аГБшт У(як + (-2 + ^ )" - ^к)
к 1 V 2/ к 2и к
к=7
Дифференцируя (25) по т и g, получим:
А ( 1 N Л 2 ( N N N Л
Ат У (як+тЬк + ^)2 )=N-7 (.У яА +тУ "к2+g у ^),
А ( 1 N \ 2 N N N Л
(. N-7 У=7 (як + тЬк + *ск )2 )=^ЫУЯл + тУ "кск + ^).
Приравнивая выражения к нулю и решая полученную систему нормальных алгебраических уравнений, будем иметь:
= У як"к У"кСк-У якСк УК
g = У Ск2 У"2 - (У"кСк )2 , (26)
ZaA Zc2-Z akck Zbk (Zbc )2-Zb2 Zc2
(27)
Аналогично случаю с моделью (2) подставим в (26) ак в явном виде из исходной модели авторегрессии и, находя математическое ожидание, выразим смещения оценок. Смещение оценки g * = X* Х2 равно
M
Zd
-gfk )bk Zbkck -Z(dk + mek + gfk )ck Zb>
Zc2 Zb2 - (ZbkCk )2
(28)
Смещение оценки m' = X* + X2 определит формула
M
Z (dk
-gfk )bk Z ck-Z (dk
-g^ )ck Zbk
(Zbkck) -ZciZb>
(29)
Как и в случае с моделью (2), величина Ьк и ск в (28) и (29) мала при относительно малом отличии друг от друга соседних пяти уровней ряда и, соответственно, велика при значительном их отличии. В отношении стохастической компоненты то же можно сказать о величинах йк, ек,/ при соседних семи, пяти и трех наблюдениях соответственно. Величину смещений (28) и (29) можно также уменьшать путем применения приема прорежения выборки. Прорежение одновременно со смещением уменьшает дисперсию оценок т' и g':
D[ g. ]=D[ Z(dk + mek + gfk- mbk- gck )bk Zbkck - Z(dk+mek+gfk- mbk- gck )ck Zbk ]
Zc2 Zbk2 -
= D
Z(dk
+ mek +
-g^ )bk Zbkck -Z(dk
+ mek +
-(Zbkck)2
-gfk )ck ZK
Z ck Zb2 - (Z bkck)2
c
c
л 85
e
45 -I----------1------------1------------1-----------1------------1-----------1
0 100 200 300 400 500 600
Номер наблюдения
---- Цена FOB на нефть Brent (дол./баррель)
.... Моделирование без прорежения (коэффициент детерминации 0,100)
----Моделирование с прорежением (коэффициент детерминации 0,792)
Рис. 2. Моделирование выборки суммой линейного тренда и двумя гармониками без прорежения
Аналогично выражается дисперсия и для m*. Оценка оставшихся параметров модели (20) предполагает подстановку в нее обеих полученных оценок частот, составление и решение СЛАУ шестого порядка для параметров линейного тренда, амплитуд и фаз гармоник аналогично первому случаю.
Эксперимент на выборке цен на нефть Brent моделью суммы линейного тренда и двух гармоник с шагом четыре прорежения и без прорежения подтвердил эффективность предложенного приема прорежения (см. рис. 2).
В результате предложенного приема прорежения выборки достигается значительное уменьшение (которое определено полученными формулами) автоковариации ошибок и £j (i Ф j), что делает возможным применение МНК, а также минимизируются смещение и дисперсия оценки частоты в моделях.
Предложенные методы позволили, в отличие от большинства известных, определять частоты (в общем случае некратные) и фазы гармоник, удобные для осуществления прогнозирования рядов.
Источники
1. Чураков Е. П. Прогнозирование эконометрических временных рядов : учеб. пособие. М. : Финансы и статистика, 2008.
2. Эконометрика / под ред. И. И. Елисеевой. М. : Финансы и статистика, 2005.
3. Семёнычев В. К., Семёнычев Е. В. Информационные системы в экономике. Эконометрическое моделирование инноваций : учеб. пособие. Самара : Изд-во Самар. гос. аэрокосм. ун-та, 2006.
4. Айвазян С. А., Мхитарян В. С. Прикладная статистика: теория вероятностей и прикладная статистика. М. : ЮНИТИ-ДАНА, 2001.
5. Poznyak A. S. Identification of non-stationary ARMA models based on matrix forgetting // Computation y Systemas. 1999. Vol. J. No. 1.
6. Сергеев А. В. Методика исследования точности обобщенных параметрических ARMA-моделей и метода их сглаживания на коротких выборках // Вестник Самар. муницип. ин-та управления. 2008. № 7.