Устойчивость в статистическом прогнозировании временных рядов Текст научной статьи по специальности «Математика»

№ 1 2006

Ю. С. Харин

Устойчивость в статистическом прогнозировании

временных рядов

1. Введение

огласно энциклопедическому словарю Уэбстера [1], «прогнозирование — это деятельность, на-

правленная на вычисление или предсказание некоторых будущих состояний или условий в ре-

зультате анализа соответствующих данных». Приведем примеры актуальных прикладных задач прогнозирования в экономике, финансах и бизнесе [2, 3]:

- разработка и анализ экономических программ роста и государственного регулирования национальной экономики на основе «больших» и «малых» эконометрических моделей;

- прогнозирование ВВП, цен и других показателей при составлении бюджета;

- прогнозирование финансовыми экспертами процентных ставок для управления денежными потоками;

- прогнозирование специалистами по маркетингу объемов продаж произведенной продукции и потребляемого сырья;

- прогнозирование «надежности заемщика» в банках;

- прогнозирование курсов акций, курсов обмена валют аналитиками на рынке ценных бумаг:

- прогнозирование состояния страховой компании.

Многие прикладные задачи прогнозирования в экономике, финансах, бизнесе и других областях сводятся к одной из важнейших проблем математической статистики — прогнозированию случайных процессов с дискретным временем иначе называемых случайными последовательностями, или временными рядами. Одним из первых исследователей, проявивших интерес к этой проблеме, является А. Н. Колмогоров [4, 5]. Математическая суть проблемы прогнозирования весьма проста: оценить будущее значение хГ4т е/т"1 в общем случае ¿-мерного случайного процесса на х е N шагов вперед по наблюдаемой в течение Т е N единиц времени реализации этого процесса

В истории исследования этой проблемы можно выделить два этапа. Первый этап, продолжавшийся до 1974 г., был ориентирован на разработку статистик (алгоритмов, или процедур) прогнозирования, минимизирующих риск (среднеквадратическую ошибку) прогнозирования для ряда фиксированных, простейших моделей временных рядов (стационарный временной ряд с фиксированной спектральной плотностью, стационарный временнбй ряд с трендом из заданного семейства функций, авторегрессионный временнбй ряд, временнбй ряд типа «скользящего среднего» и др.). Обзор этих результатов имеется в [6-8].

В 1970-е годы обнаружилось, что многие из разработанных «оптимальных» процедур прогнозирования при решении прикладных задач имели риск, значительно больший ожидаемого теоретического значения. В 1974 году швейцарский математик Питер Хьюбер в докладе на Математическом Конгрессе в Ванкувере [9] сформулировал причину этого явления: «Статистические выводы (в том числе статистические прогнозы) лишь отчасти базируются на наблюдениях. Такой же важной их основой являются априорные предположения, совокупность которых образует гипотетическую модель исследуемого процесса МП». В прикладных задачах гипотетические модельные предположения М0 часто нарушаются и это может приводить к неустойчивости прогнозирующих статистик, оптимальных только для модели М0. Известны следующие основные типы искажений гипотетической модели М0 [10-13]: «негауссовость» ошибок наблюдения; зависимость

{*1.....хг}сЯ".

82/

Ne 1 2006

(или коррелированность) ошибок наблюдения; неоднородность ошибок наблюдения; наличие | «выбросов»; ошибки спецификации модели; наличие «разладок»; наличие «пропусков» в наблю-даемом временном ряде. Было предложено [10] строить так называемые устойчивые, или робас- <j тные (от англ. robust — крепкий, стойкий, устойчивый) статистические выводы — это статисти- g ческие выводы, «слабочувствительные к малым искажениям гипотетической модели М0». С этого момента начинается второй этап исследования проблемы статистического прогнозирования временных рядов. В странах СНГ развитие робастной статистики было поддержано научной школой С. А. Айвазяна [2, 3].

Отметим, что большинство публикаций по проблеме робастности в статистическом прогнозировании посвящены в основном задачам статистического оценивания параметров модели и статистической проверки гипотез при наличии искажений. Хотя эти задачи — необходимые элементы проблемы робастности в статистическом прогнозировании, они не покрывают исследуемую проблему полностью. Актуальны следующие задачи, которым и посвящается данная статья:

а) математическое описание искажений гипотетической модели М0;

б) количественная оценка устойчивости традиционных прогнозирующих статистик при наличии искажений;

в) построение новых робастных статистик прогнозирования.

2. Вероятностная модель и характеристики устойчивости прогнозирования

Пусть наблюдается d-мерный временнбй ряд хрК1 с дискретным временем teZ,X={x^...,xT'),sRTd— составной вектор наблюдений для Г моментов времени (штрих означает транспонирование), TeN — длительность наблюдения временного ряда, xT+tsRd — ненаблюдаемый случайный вектор, подлежащий прогнозированию в будущий момент времени Т+х, xeN; в эконометрических приложениях Т называется базой прогнозирования, at — горизонтом прогнозирования. Вероятностная модель наблюдаемого временнбго ряда при наличии искажений определяется семейством вероятностных мер

{Р\еЛА),А е BTd : Т е А/, 6° е 0 с /Г,е е [0,е+]},

где BTd —борелевская о-алгебра подмножестве RTd; 9° — неизвестное истинное значение параметров модели временнбго ряда; е — уровень искажений модели, а е+ > 0 — его максимально допустимое значение. Если е+ = 0, то искажения отсутствуют, и мы имеем гипотетическую модель MQ. Статистический прогноз xT+t е Rd определим с помощью статистики

а его погрешность — с помощью среднеквадратического риска (среднеквадратической ошибки) прогнозирования

pE=pe(^;9°) = E{llxr+t-x7+Tll2}>0, (1)

где Е{-} — символ усреднения по мере у С^Д). Для гипотетической модели М0 (когда е = 0) функционал (1) будем называть гипотетическим риском: р0 = PoCrv^0)-

Чтобы устранить зависимость от неизвестного параметра 9° в (1), мы будем использовать также интегральный риск:

r£ = ^ = Jpe(k;e)it(e)c/8;>0, (2)

0

где л(9) — некоторая подходящая весовая функция, например, плотность распределения вероятностей, если 9° — случайный вектор. Заметим, что для некоторых моделей временных рядов, например,

83

Ns 1 2006

i-

i ss

u to

¡2 o

o §

для регрессионных моделей (см. разд. 3), точечный риск (1) не зависит от 0°, так что интегральный и точечный риски совпадают: лЕ = рЕ.

Гарантированным (верхним) риском называется точная верхняя грань множества значений риска при всевозможных допустимых искажениях модели М0:

г+ = r^fu) = sup0<E<£i ге(У ■ (3)

Пусть x°+t = /£(X;e°) — оптимальная прогнозирующая статистика при известном 9° и отсутствии искажений (е = 0), которая доставляет минимальное значение р0 = p0(f°vQ°) точечного риска Ро = Ро(^;е°) и интегрального риска г0 = |еро(/7°;0)л(0)с/0, причем л0 >0 за исключением вырожденных случаев. Так как 0° неизвестно, то для прогнозирования используется семейство так называемых подстановочных («plug-in») прогнозирующих статистик:

(4)

где бе Rm — некоторая состоятельная оценка параметра 0°. Если limr_>„(po(/r-T;0o)-po) = O, то прогнозирующая статистика (4) называется асимптотически оптимальной. Следуя [12,14], относительное приращение гарантированного риска

к = к(у = (/-+(у-д/л0>0 (5)

назовем коэффициентом неустойчивости риска, а величину

е" = е'(§) = suple: к(у < 5} — (6)

8-допустимым (5>0) уровнем искажений. Величина е* представляет собой максимальный уровень искажений, при котором коэффициент неустойчивости риска не превосходит заранее заданного значения 8. Чем меньше к и больше е", тем более робастной (устойчивой) является статистика прогнозирования fTz(■) по отношению к рассматриваемым искажениям. Поэтому прогнозирующую статистику f^ ^ назовем минимаксно устойчивой, если минимален коэффициент

к^Нпи.к (а (7)

§ !

3 неустойчивости:

I

I

в В [11] предложена еще одна характеристика «качественной робастности» — «пороговая точка» е

| Хампеля — это максимально допустимая доля «выбросов» в выборке, при произвольном варьирова-

§ нии значений которых прогнозирующая статистика остается ограниченной.

I

§

х

3. Устойчивость прогнозирования регрессионных временных рядов

о 3.1. Случай функциональных искажений

с"

о Пусть наблюдаемый ¿/-мерный временной ряд х, е Я" определяется стохастическим уравнением: |

где ¡;( еЯй — случайный вектор ошибок наблюдения в момент времени ц е(7сЯм — заданный вектор М «независимых переменных» («регрессоров») в момент с,

и - множество допустимых значений регрессоров, {\|/.(«):/?м-»/?} — некоторая система т линейно независимых функций, такая что det(Zгt=lv|/(uДl|/(utЮ*Ol \|/(1/)=(\|л(и))е/Г, е°=(0о(у)=(е°(]):....:.ео(т))е^><п1-(с^хт) — | матрица неизвестных истинных значений коэффициентов регрессии,

2 6°, еЯ" — /'-й столбец матрицы е°Д(-):Ям ->ЯС' — некоторая неизвестная неслучайная векторная функция, описывающая многомерные функциональные искажения (ФИ) гипотетической модели. Слу-

84

Нв 1 2006

чайные ошибки {£,} предполагаются независимыми в совокупности, одинаково распределенными, | причем ^

Е£,} = 0, Сс«/Й„$,} = 1, с1е12:*0. Ч

г

При Х(и,) = 0 (8) превращается в гипотетическую модель М0 — многомерную (множественную при с/ = 1) линейную регрессионную модель временного ряда [7] с гипотетической функцией регрессии Ци) = у01(и)) = в\(и). Если М = \и, ^¿,то получаем трендовую модель временного ряда [6-8], наиболее часто используемую при прогнозировании (см. ниже пример 1).

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

ФИ-1. Интервальные искажения:

£,» < Х,(и) < £1+(и),и е и,1 = 1.....б,

где е±(и) = (е/±(и)) — некоторые заданные граничные векторные функции. ФИ-2. Относительные искажения:

I Х,-(и) I /1 1< е,,е, > 0,/' = 1.....с/.иеЦ

где заданный вектор е = (е,) определяет относительные уровни искажений для компонент многомерной линейной функции регрессии; относительная «ошибка спецификации» /-Й компоненты не превосходит е, • 100%.

ФИ-3. Искажения в /, -метрике:

I11 II, 1ЦиТ„)Ые,

где е>0 —заданный уровень искажений.

Исследуем вначале устойчивость традиционно используемого алгоритма прогнозирования временного ряда (8) при искажениях ФИ-1, ФИ-2, ФИ-3 с помощью характеристик устойчивости (1М7).

Примем обозначения:

А = I)'V'(Ц)е Ятхт, С = Хх,\|/'(и,)е Я*™,

(=1 (=1

I- = ея*™, ег+Т=м-уое яа,

(=1

а, = у'{и,)АМит+Л * = 1.....Т, К0 = <(и7+Т)Л->(иг+Т) > 0,

(г)+ = тах(г.О):

Прежде всего отметим, что при отсутствии искажений {Х(и,) = 0) и известном 0° минимальный риск

л0 = р0 = ?л(Х)

обеспечивает оптимальная прогнозирующая статистика

При решении прикладных задач прогнозирования с неизвестным 0° традиционно используется предположение IАI* 0 и алгоритм МНК (метода наименьших квадратов), основанный на подстановочной прогнозирующей статистике (4):

Хг+х=0у(О,0 = СЛ-1, (9)

где 0 = (0(1) :...:0(т,) 6 Я"хт — МНК-оценка матрицы 0°.

85

N9 12006 ЧЧчЧ

ПРИКЛАДНАЯ ЭКОНОМЕТРИКА

Прямым вычислением по формулам (1), (2), (8), (9) с использованием введенных обозначений легко показать, что если наблюдаемый временнбй ряд удовлетворяет многомерной линейной регрессионной модели с искажениями (8), то риск прогнозирования (1), (2) для МНК-алгоритма (9) имеет вид:

ле(7д) = ре(Гд) = (1 + К0)г0+ И - Чит„) II2.

Отметим, что если «план эксперимента» (ц.....ит}аи выбран так, что минимальное собственное значение при 7"-»°°, прогнозирующая статистика (4) является асимптотически оптимальной:

/<0-> 0, л0(7"д) = (1 + К0)л0 -> г0.

Теорема 1. Если имеют место искажения ФИ-1, то гарантированный риск (3) для МНК-алгоритма (9) равен

;(7",т) = г0(Т,т) + 1тах(хЦ(а,1е,,±(и,)- (-а,)+е^(и()) - е,т(иг+х))2

Следствие 1. Коэффициент неустойчивости риска (5) равен

к(7д) = К0 + ^ £ тах (х[=1 ((а,)+ ек± (и,) - (-а,)+ект (и,)) - е^ (иг+х ))2 ^ /(?/"Е).

Следствие 2. В случае множественной линейной регрессионной модели (с/ = 1,0{£(} = а2) с искажениями ФИ-1 коэффициент неустойчивости риска равен

к(Т,%) = К0 + ^г(иЪа, ГГ.

а V м )

Из этого выражения следует, что коэффициент неустойчивости риска содержит два слагаемых. Первое слагаемое К0 ->0 при и не зависит от искажений ФИ-1; оно порождено статисти-

^ ческой погрешностью состоятельного оценивания коэффициентов регрессии 9° по конечной вы-борке объема Т без учета ФИ-1. Второе слагаемое порождено искажениями, оно тем больше, чем | больше вклад систематической ошибки е+2 в уравнении (8) по сравнению с дисперсией случайной 5 ошибки а2.

| Аналогично Теореме 1, с учетом ограничений на (Х(и()}, задаваемых ФИ-2, можно установить, что « если имеют место искажения ФИ-2, то гарантированный риск (3) для МНК-алгоритма (9) равен

| г+(Гд) = г0(Гд) + ё'х,

9

I

£ ё = (е2) 6 Яб,х = (Х„) 6 Я'.Х* = (sr.il (еЧ(ч))* II а, I +1 (9>(иг+х))к |).

§

При этом коэффициент неустойчивости риска

| к (Гд) = К0 + ё'х/(т §

§ если к тому же уровни искажений одинаковы (е, = ... = га =е), то 8-допустимый уровень искаже-

| ний

Е е+(6) = ((6-/<0)+(ГлХ)/1хкЧ1/г

0 V

19

¡2 Если имеют место искажения ФИ-3, гарантированный риск МНК-алгоритма (9) равен [15]:

1 г+ (Г д) = г0(Г д) + е2 тах{1, а2,..., а2).

I

Л)

При этом справедливы следующие выражения для коэффициента неустойчивости риска и 8-до-

пустимого уровня искажений (8 > 0):

Ив 1 2006

к(7",т) = К0 + е2тах{1,а12.....а2}/(?л2> |

е+(8) = ((*лГ)(8 - /<0)+)1/2/тах{1,а2.....а2 }• ¡5

О

Теперь рассмотрим задачу построения робастных прогнозирующих статистик для случая интер- д вальных симметричных искажений ФИ-1: е_(и) = -е+(и) Будем использовать подстановочную прогнозирующую статистику:

хт„ = Щит„),

отличающуюся от (9) использованием оценки из семейства так называемых /М-оценок [10, 11] вместо МНК-оценки 9;

0 = агдггнпе1р(х(-0у(и,)), (10)

(=1

где : Rd —»[0, +«>) — некоторая выпуклая четная относительно 2 функция потерь, монотонно возрастающая относительно I г, 1,...,1 zd I.

В случае искажений ФИ-1 нами предложен следующий выбор этой функции потерь в (10):

p(z) = Ip1(z/), Pl(z,):

/=1

(z,+ßei+(u,))2, z, < —е,, 0, lz,.l<£

(z,-ße;+(iv,))2, z, > e,,

где ß > 0 — некоторый параметр.

Заметим, что функция потерь такого вида при ß = 1 получается, если ошибки наблюдения {£,} предполагаются гауссовскими и используется метод максимального правдоподобия относительно

расширенного набора параметров 8,{X(u,):f = 1.....7"} при ограничениях, заданных ФИ-1. Методом

Монте-Карло строились оценки коэффициента ß из условия минимума коэффициента неустойчивости риска [15]:

к min-

и

Как показали эти эксперименты, оптимальное значение параметра ß* = 1.

Отметим также, что вычисление оценки 0 сводится к задаче выпуклого программирования, которая решается с помощью стандартной итерационной процедуры. Для увеличения скорости сходимости этой процедуры функция потерь p(z) сглаживается с помощью кубической параболы [15].

3.2. Случай аддитивных «выбросов»

Пусть теперь наблюдается одномерный (d = 1 для простоты выкладок) регрессионный временной ряд с аддитивными искажениями-«выбросами» [10,11,13,15]:

х,=(в°)>(и,) + iiA+^fe/V, (11)

где 0° = (0°.....Q°m)'eRm — вектор-столбец неизвестных коэффициентов множественной линейной

регрессии;

т|, е{0,1} — случайная величина Бернулли, Р{г), = 1} = 1-Р{т|( = 0} = е;ее[0,е+] — вероятность появления «выброса», 0<е+ <1/2;{й,} — независимые одинаково распределенные случайные величины, {fy} = a,D{i},} = Кс2,К > 0;{^(},{т|,},{й(} — независимы в совокупности.

Обозначим: lm —единичная (тхт) — матрица;

0Т,\ е RT — векторы-столбцы, все компоненты которых равны 0 и 1 соответственно; Уг = (у;(ц)) е RTxm;XT = (х„...,х,)' е RT;g(T, т) = ^(Ч^Ж J e R7-

87

Не 1 2006

МНК-статистика прогнозирования на х шагов (построенная для гипотетической модели при е+ = 0) является частным случаем (9) для d = 1:

х7+т = é>(ur+t),e = OFV 4V)-14"r xr (12)

Проводя эквивалентные преобразования с помощью (1)-(3), (11), можно получить следующее утверждение.

Теорема 2. Пусть наблюдаемый временной ряд {х,} удовлетворяет регрессионной модели (11) с «выбросами», причем IA l=lf'rTT |*0. Бели для прогнозирования случайной величины xr+t используется МНК-алгоритм (12), гарантированный риск (3) равен

л+(Г,т) = (о2 + е+(а2 + Ко2))( 1 + К0) + е2а2((1 - Гт д(Т,х))2 -К0~ 1).

Следствие 1. В случае «выбросов по дисперсии» (а = 0,К>0) коэффициент неустойчивости (5) и 8-уровень искажений имеют вид:

к(7-,т) = К0 + е+К(1 + К0),

e+(5) = m¡n{1/2,(5-/<0)t(/<(1 + /<0))-1}.

Следствие 2. В случае «выбросов по среднему» (К = 0,а * 0) коэффициент неустойчивости (5) равен

к(7",т) = К0 + е+(а/о)2(1 + К0) + е!(а/о)2(( 1 -1'т <?(Г,т))2 - К0 -1).

Для повышения устойчивости регрессионного прогнозирования при наличии «выбросов» целесообразно использовать локально-медианный подход [16]. Пусть 3 = {1,...,7"} — множество моментов времени, 3' = {?{",...,t(n"} с 3 — 1-е подмножество п упорядоченных в порядке возрастания моментов времени (/ = 1.....L),L — число рассматриваемых различных подмножеств (1</_<С"), п — число моментов в каждом из подмножеств 3' (т<п<Т). Определим соответствующие подматрицы матриц e Хт,*¥т : X''1 = (х((„) е R", Ч^1 = (v|/у(и„,)) е Rnxm, а с их помощью аналогично (12) — 1-ю локальную МНК-

О \ ' / i ,

| оценку 9(,) = ((Т(п")'Ч/''))-1(У^),Х'') и 1-й локальный прогноз Xj^ = (é(") \|/(uT+t). Локально-медианным

5 (ЯМ) прогнозом называется выборочная медиана локальных прогнозов [16]:

| XT+%=med^L (13)

Отметим, что если п = Т, то ЛМ-прогноз (13) совпадает с МНК-прогнозом (12).

| Доказано [16], что если = С", имеет место модель с «выбросами» (11) и используется ЛМ-прогноз §

о. уравнения:

(13), «пороговая точка» Хампеля е е [0,1-л/7"], и на этом промежутке является единственным корнем

V, t П0-7)-

(=0 /

8. т У >- ,

с При этом, если размер подвыборки п фиксирован, а Т -> о»,то е 1-2 1/п, и максимум устойчивости (в смысле «пороговой точки» Хампеля) достигается при оптимальном размере подвыборок:л" =т. Таким образом, размер п каждой локальной подвыборки {1 = \...,0 целесообразно задавать равным числу оцениваемых параметров т.

Асимптотические разложения гарантированного риска и коэффициента неустойчивости риска получены в [16].

Пример 1. Для сравнения традиционного МНК-прогноза (12) и ЛМ-прогноза (13) рассматривалась прикладная задача прогнозирования динамики объема выручки от реализации продукции Новопо-| лоцкого производственного объединения «Полимер» концерна «Белнефтехим» на основе частного случая модели (1):

1

£ d = 1, m = 3, \|/(f) = (1,f,í2)', 7=14, те{1,2,3};

88

Не 1 2006

при этом х, — это логарифм объема выручки. Использовались реальные ежемесячные статистичес- |

кие данные за 14 месяцев, по которым необходимо было построить прогноз на последующие 1-3 ^

месяца. Графики исходного временного ряда, МНК-прогноза и ЛМ-прогноза (с параметром п = т = 3 ^

) представлены на рис. 1. Модуль ошибки прогнозирования объема выручки для ЛМ-метода (13) по д сравнению с МНК-методом (12) оказывается меньше приблизительно в 4 раза.

4. Робастность прогнозирования авторегрессионных временных рядов

4.1. Случай функциональных искажений

Пусть наблюдаемый временной ряд х, еЯ1 определяется стохастическим разностным уравнением т-го порядка:

х,=(в°уХ,_,+Ц0 + ^, teN, (14)

где в° = (в°)ейт —вектор-столбец т авторегрессионных коэффициентов; Х,_,=(х,_,.....х,_т)'еНт,

Х0=От, — независимые в совокупности гауссовские случайные величины с распределением Л/(0,<т2), %{■) — некоторая неизвестная неслучайная действительная функция, описывающая функциональные искажения гипотетической модели АР (т) (гипотетическая модель АР (т) получается из (14) при Я.(-) = 0).

Как известно [14], оптимальный прогноз для гипотетической модели АР (т) при известном 9° определяется рекуррентным соотношением:

x0T+T = (Q°yx°T+hVi = ^.....Т, 05)

где Х°т+И = (Х°т+Н.....Х%_ту- = Х5 для 5 < Т.

Обозначим (тхт) — блочные матрицы:

/ (9°)' 4 ( 9' >

В0 = у/т-1 0^.1 , ,в = к 'т-1

89

Ив 1 2006

б,,(б,) — /'-я строка (/'-й столбец) матрицы В, (б)и — (/с,/ )-й элемент матрицы В. Будем предполагать, что выполнено условие устойчивости модели АР (т): все собственные значения б0 лежат внутри единичного круга.

Используя [7] и (1), (14), (15), можно показать [16], что если наблюдаемый временной ряд х( соответствует гипотетической АИ (т) — модели (14) при наличии функциональных искажений >.(•) и присутствует «ошибка спецификации модели» (т.е. вместо истинного значения 6° в алгоритме прогнозирования (15) используется значение 9 * 9°), в этом случае риск прогнозирования (1) равен

р0(9°,в;Т,х) = о^1 + | ((В^Л^-ВЖХ)2)- 06)

Если присутствует «ошибка спецификации» 0-0° е/?т и функциональные искажения ограничены:

Т Т+х

I 07)

(=1 (=.741

где е(1),е(2) >0 — уровни искажений для интервала наблюдения и интервала прогнозирования соответственно, тогда гарантированный риск прогнозирования

р+(9°,9;7",т) = р0(9°,9;Т,х) +

+| £(1)^7~Х(((Вг -В1)В'0)и)2 +е(2)утХ((бо)11)2

Рассмотрим теперь ситуацию, когда в (15) вместо неизвестного истинного значения вектора коэффициентов авторегрессии 9° подставляется случайный вектор § = (9.) е дт —статистическая оценка 9° по методу наименьших квадратов [7]:

| ё = С I I (18)

| построенная на основе временного ряда X = (х1,...,х7-о)еЯт°, наблюдаемого (на «этапе обучения») не-

5 зависимо от временного ряда X = (хг...,хТв)е Ит.

^ Пусть прогнозированию на один шаг (х = 1) подлежит временной ряд (14), соответствующий мо-

6 дели АИ (т) с искажениями, удовлетворяющими (17) при е(1) = е+ > 0, е(2) = 0, и для прогнозирования ис-| пользуется «подстановочная» статистика (15):

О Хт+1=в'Хт,

8 где оценка 9 определена (18) и имеет матрицу вариаций Е{(9-9°)(9-90)'} = (7"оГ1Р, Р —некоторая ^ симметричная положительно определенная (тхт) -матрица. Тогда гарантированный риск прогнози-рования равен [16]

| р+(9° Т,7-0) = ог + ^1(бГ)'1Р(0, + е^цтах(Сг),

^ 'о 'о

| где ^(б,-) — максимальное собственное значение (Гх7")-матрицы Сг=^(бо~1) ^(б^1).^, /,у=1,...77

0 4.2. Прогнозирование векторных авторегрессионных временных рядов ¡а при наличии пропущенных значений

§ Пусть определен ¿/-мерный (с/>1) временной ряд, удовлетворяющий гипотетической модели

1 \/АЯ(1) векторной авторегрессии первого порядка [7,13,17]:

| У, = бУм+Е„ 09)

На 1 2006

где У, = (у(1,...,уи)'еЯ" —случайный вектор наблюдения в момент времени инновационный про- | цесс {е, = е Я" : ? € 1\ — независимые одинаково распределенные гауссовские случайные

векторы, ^

Е{Н,} = 0,, Е{5,Н'(} = 1, IЕ I* 8

б = (Ь,) — (с/х ^-матрица коэффициентов авторегрессии, все характеристические числа которой лежат внутри единичного круга. Заметим, что рассмотренная выше модель АН(.т), а также векторная авторегрессия \/АЯ (т) сводятся к исследуемой здесь модели Х/АИ (1) расширением пространства состояний Я" [7,13].

Искажения модели (19) заключаются в наличии пропущенных значений при наблюдении процесса {У,}. Для описания таких искажений введем в рассмотрение детерминированную векторную функцию О, = (оп,...,ош), ? е которую назовем шаблоном наблюдений:

[1, если значение у,, наблюдается, о,, = <!

[О, если значение уе пропущено.

Примем обозначения: М = {(?,/'),? |1.....с1}: о„. =1} — множество лексикографически упорядоченных пар индексов наблюдаемых значений временного ряда (19): К=\М\ — общее число наблюдаемых значений;^ = плю{?:£?=1ой >0), Т = тахЦ :£?=1о(; >0}— начальный и конечный моменты времени регистрации наблюдений; к = = %(к) — биекция: М <-> {1.....К}; X = (х,.....хк)' е Як — составной вектор

всех наблюдаемых значений; УТ+Х е. Я" — подлежащий прогнозированию вектор в будущий момент времени Т + %; Р = Соу{Х,Х}еЯКхК, в = Соу{Х,УГи( е ЯКх" , /-/= Соу{У7+т,У;+т} е Я"*0' — ковариационные матрицы.

Известно [7], что в случае полных данных (К = б(Т - +1)) и известных параметров 6,Е модели (19) минимальный риск р^ = ¿/"(Х^В' 1(6')') > 0 достигается для прогнозирующей статистики У°+т = Б%

Имеют место следующие выражения для ковариационных матриц Р = Щ), в = (д^), Н = (Л. )

= 1.....к,1>г,

д1/={вт-^н)шп,1 = 1.....к,у = 1.....а-

н = 1БТ(е'у.

/=о

Теорема 3. Если имеет место модель (19) с «пропусками», параметры Б,Х известны и 1Р1*0, то прогноз по критерию максимума правдоподобия и его риск определяются соотношениями:

Ут» = I = Л>(ВЯ)ХЛ(В,1) = (20)

р 0=ПН)-НвТ-'в). (21)

Следствие 1. Коэффициент неустойчивости риска, порожденной «пропусками», равен

к0 = = Г£/Л(в' 1(В'У)-К(вГ-'в))/Ъг (В1 Х(В'У) > 0.

р0 Ч'=Х ) /=о

Следствие 2. Если имеет место модель (19) с «пропусками», IРI* 0 и параметры неизвестны, то прогноз по критерию максимума правдоподобия является «подстановочным»:

V'т„=А0ф,±)Х. (22)

где В,1 — оценки максимального правдоподобия, являющиеся решением экстремальной задачи:

ЦВЛ) = Х'Р"1Х + 1п IРI +1п IН -вТ-'в 1-> гтш- (23)

е,1

Из-за вычислительной сложности задачи (23) предлагается в (22) вместо оценок максимального правдоподобия использовать МНК-оценки Андерсона [7], которые модифицируем применительно к ситуации с «пропусками» следующим образом:

Не 1 2006

где 6 и Я вычисляются с помощью статистик:

.....*

Некоторые результаты исследования предложенных оценок и риска прогнозирования представлены в [19, 20].

Пример 2. Для численного анализа влияния «пропусков» на риск прогнозирования рассматривалась ставшая классической в зарубежной учебной статистической литературе задача статистического прогнозирования размера популяции канадской рыси [18]. На основе реального (центрированного) временного ряда ежегодных размеров популяции рыси, собранного в течение 7 = 113 лет, построена модель авторегрессии порядка с1 = 11:

у, =1,0938ум -0,3571у,_2 -0,1265у,_4 +0,3244умо -0,3622уМ1 +$„

= о2 =0,04405.

со

0

! ¡5

1 I

Я

I

г

§

и ¥ 5

и $

и в

Й и

0

1

>5

8 £

В этот временной ряд искусственно внедрялся единственный «пропуск» в момент времени те{7-с/ + 1,...,7-1}, и рассматривалась задача прогнозирования значения уг+1 на основе (у,: ? = 1.....т-Хт + 1.....7} при двух уровнях априорной информации:

а) коэффициенты авторегрессии известны;

б) эти коэффициенты неизвестны.

Вышеприведенная модель АР (11) преобразуется в модель (19) для векторного временнбго ряда

П=(У„Ум.....Умо)'еЯ11.

На рис. 2 для случая а) изображен график теоретической зависимости коэффициента неустойчивости риска к(т), вычисленного согласно (20), (21) и следствию 1, а также приведены 95%-е доверительные интервалы для к (т), вычисленные методом Монте-Карло. В случае б) границы доверительных интервалов меняются лишь на 1-2% (доля «пропусков» е = 0.01 достаточно мала). Из рис. 2 видно, что местоположение «пропуска» во временном ряде оказывает существенное влияние на величину риска прогнозирования. Так, например, «пропуски» в моменты времени 7-1,7-9,7-10 более значимы, чем в другие моменты времени.

5. Заключение

Результаты, полученные в данной обобщающей статье,

$ теор.

1 знач.

1 верхи.

I —в— довер.

1 гран.

1 нижн.

а! —й—довер.

Лл гран.

Т-10 Т-* Т-4 Т-2

Рис. 2. Коэффициент неустойчивости риска при наличии «пропусков»

Не 1 2006

предоставляют исследователю количественные оценки гарантированного среднеквадратического | риска прогнозирования (3), коэффициента неустойчивости риска (5) и 5-допустимого уровня искаже- ^ ний (6) для существующих методов статистического прогнозирования и типовых искажений («ошиб- ^ ки спецификации», функциональные искажения, «выбросы», «пропуски») трендовых, регрессионных и д авторегрессионных моделей наблюдаемых временных рядов. Эти оценки позволяют раскрыть механизм влияния различных искажений на риск прогнозирования и указать пределы изменения уровней искажений, при которых обеспечивается заданная надежность прогнозирования. По минимаксному критерию (7) построены робастные прогнозирующие статистики и найдены условия их эффективного применения. Методы асимптотического анализа традиционных прогнозирующих статистик и синтеза робастных статистик могут быть распространены на другие модели временных рядов и другие типы искажений.

Представленные в статье теоретические результаты нашли подтверждение в численных экспериментах, проведенных как на модельных, так и на реальных данных, и были внедрены в пакеты прикладных программ РОСТАН (РОбастный Статистический АНализ), СТАТПРО (СТАТистическое Прогнозирование) и СЭМП (Система Эконометрического Моделирования и Прогнозирования), разработанные в Белорусском государственном университете.

Исследования были частично поддержаны грантом ИЧТАБ № 03-51-3714.

1. Maclaren A. Consolidated-Webster encyclopedic dictionary. N.-Y.: СВР, 1946.

2. Айвазян С. А. Основы эконометрики. М.: ЮНИТИ, 2001.

3. Айвазян С. А. и др. Прикладная статистика, т.1-3. М.: Финансы и статистика, 1983,1985,1989.

4. Колмогоров А. Н. К вопросу о пригодности найденных статистическим путем формул прогноза/ /Журнал геофизика. 1933.1 3.

5. Колмогоров А. Н. Интерполирование и экстраполирование стационарных случайных последова-тельностей//Изв. АН СССР. Серия материалов. 1941. Т. 5.

6. БоксДж„ Дженкинс Г. Анализ временных рядов. Прогноз и управление. Вып. 1, 2. М.: Мир, 1974.

7. Андерсон Т. Статистический анализ временных рядов. М.: Мир, 1976.

8. Bowerman В. L, O'Connel R. Т. Forecasting and Time Series. Belmont: Duxbury Press, 1993.

9. Huber P. Some Mathematical Problems Arising in Robust Statistics//Proc. Int. Congress of Math. Vancouver. 1974.

10. ХьюберП. Робастность в статистике. М.: Мир, 1984.

11. Хампель Ф. Робастность в статистике. М.: Мир, 1989.

12. Kharin Yu. Robustness in Statistical Pattern Recognition. Dordrecht: Kluwer, 1996.

13. Greene W. H. Econometric analysis. N.-Y.: Macmilllan, 2000.

14. Kharin Yu. Robust Forecasting of Parametric Trends//Studies in Classification and Data Analysis. 2000.

15. Харин Ю.С., Маевский B.B. Робастность регрессионного прогнозирования при наличии функциональных искажений модели//Автоматика и телемеханика. 2002. № 11.

16. Kharin Yu. Robustness Analysis in Forecasting of Time Series//Development in Robust Statistics. N.-Y.: Springer, 2000.

17. Abraham В., LedolterJ. Statistical Methods For Forecasting. N.-Y.: Wiley, 1989.

18. TongH. Non-linear Time Series. Oxford: Clarendon Press, 1999.

19. ХаринЮ. С., [урин А. С. Прогнозирование векторных авторегрессионных временных рядов с про-пусками//7еория вероятностей, математическая статистика и их приложения. Минск: БГУ. 2004.

20. Kharin Yu., Нигуп А. «Plug-in» Statistical Forecasting of Vector Autoregressive Time Series With Missing Values//Austrian Journal of Statistics. 2005. Vol. 34. № 2.

Литература

Vol. 17.

93