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

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2011 Управление, вычислительная техника и информатика № 1(14)

УДК 512.234

В.А. Симахин АДАПТИВНЫЕ РОБАСТНЫЕ НЕПАРАМЕТРИЧЕСКИЕ АЛГОРИТМЫ ПРОГНОЗА

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

Ключевые слова: временной ряд, робастность, непараметрический, прогноз.

Статистическому анализу временных рядов посвящено большое количество исследовательских работ (см., например, [1 - 3, 5]). В качестве математической модели случайного процесса используются общие гетероскедастические модели нелинейной авторегрессии

X = т(Хг_т-1) + g(Хг_т_1)в,., (1)

где Х,_т_1 = (X,■-1,..., X'-т )Т , т(Х,_т_1 ) - фуНКция аBTорегрессии, g(Х,_т_1) -

функция волатильности, е, - последовательность независимых и одинаково распределенных случайных величин. При непараметрической постановке задачи структура функций т(Х,_т_1), g(X,_т_1) и распределение е, считаются неизвестными, удовлетворяющими достаточно общим предположениям [2, 3]. Необходимо отметить, что работы на основе непараметрического подхода занимают достаточно скромное место в исследованиях временных рядов, тому свидетельство обзоры [2, 3]. Большая часть этих работ связана с анализом условных функционалов в условиях зависимой выборки [2, 3, 5]. Начиная с 60-х гг. прошлого столетия возрос интерес к исследованиям, связанным с созданием робастных алгоритмов прогноза временных рядов [4]. Для задачи прогноза в непараметрической постановке задачи предложено три непараметрических алгоритма - регрессионного, медианного и модального типов [2, 3, 5 - 7]. Данные алгоритмы были получены по аналогу регрессионных задач путем переноса их на авторегрессионные задачи. Алгоритм регрессионного типа [3, 5, 6] не является робастным и, очевидно, считалось, что алгоритмы медианного и модального типов решают задачу робастного непараметрического прогноза. Причем робастные свойства алгоритма модального типа мало изучены. Анализ показывает, что алгоритмы прогноза, наследуя свойства алгоритмов регрессионных задач, в то же время не в полной мере учитывают специфику задачи прогноза. Прогноз проводится по последним членам временного ряда, и любой выброс в этих данных приводит к катастрофическим результатам. Меры борьбы с выбросами в независимых переменных, применяемые в регрессионном анализе, в данном случае неприменимы. Необходимо принимать дополнительные меры, чтобы получить робастные алгоритмы прогноза. В данной работе на основе взвешенного метода максимального правдоподобия [9] синтезирован класс робастных (устойчивых) непараметрических алгоритмов для задачи прогноза. Данные алгоритмы включают как частный случай непара-

метрические алгоритмы регрессионного, медианного и модального типов, введенные ранее [2, 3, 5-7].

Рассмотрим задачу прогнозирования случайного стационарного в узком смысле процесса X(/) вида (1), удовлетворяющего условию слабой зависимости (сильного перемешивания (СП) или равномерно сильного перемешивания (РСП)) [5] с неизвестной функцией распределения Е(х,t) еП, где О - непараметрический класс абсолютно непрерывных функций распределения.

Пусть х1,...,хм - выборка из случайного процесса X(/), наблюдаемого на интервале [0, Т] (х, = X(/,), ^+1 _ ti = At, хм = X(Т)). Требуется по т последним значениям временного ряда хм_т = (хм,..., хм_т_1)Т найти прогноз XN+k на к тактов вперед.

В общем случае оптимальный прогноз в зависимости от критерия качества есть функционал от условной функции распределения Е ^ / XN _т)

Вид функционала 3 зависит от используемого критерия качества. В частности, для квадратичной функции потерь оптимальным прогнозом является условное среднее. В дальнейшем были введены алгоритмы на условной медиане Е(XN+к / XN_т) = 0,5 и условной моде XN+к = а^шах Е'(t/ XN_т) [2, 3, 7]. Дан-

ные алгоритмы рекомендовались к применению в случае, если выборка «засорена» (имеются, например, единичные дельтообразные выбросы).

В непараметрическом случае, когда вид Е(х, t) неизвестен, в качестве алгоритмов прогноза берутся оценки

где ^ / % _ т) - непараметрическая оценка условной функции распределения.

В качестве Ем^ ^ / хм^ _т) будем использовать две классические оценки из класса дискретных и непрерывных оценок, введенных в [10] и [11] соответственно:

1. Постановка задачи и анализ выбросов

1.1. Постановка задачи

(2)

ХЫ+к = ((^1 -т )) ,

(3)

N - т-к

(4)

где

(6)

оценка плотности Розенблатта - Парзена;

xi = (xixi+m-1)T , 0 < G(u) < 1, G(-да) = 0, G(^) = 1, G(-u) = 1 -G(u),

Iu2dG(u) <да, G'(u) = K(u), IK(u)du = 1, |uK(u)du = 0, |u2K(u)du <да ,

IK2 (u)du = Q < да, hN ^ 0, Ж • hm ^ да , при N ^ да ,

C(u) = [1, u > 0; 0, u < 0] - функция Хевисайда.

Замечание 1. Выше приведены классические ограничения на ядерные функции K(u), G(u) и коэффициент размытости hN . На коэффициенты перемешивания процесса обычно накладываются условия, обеспечивающие сходимость оценок,

да

например для РСП ^ VP(T) <да [5, разд. 3.1]. В дальнейшем будем считать, что

Т=1

данные условия регулярности выполняются [2, 3, 5 - 7].

Замечание 2. К настоящему времени предложены различные алгоритмы на основе активных ядер (см., например, [2, 3, 5]), позволяющие повысить скорость сходимости оценок. Применение методов получения активных ядер (например, метод локальной аппроксимации) усложняет математику, но не меняет сути, в связи с этим в данной статье на них останавливаться не будем.

1.2. Анализ выбросов

Без потери общности будем считать к=1. Основные алгоритмы прогноза относятся к «гусеничным алгоритмам». По временному ряду Xj,...,xN прокатывается

гусеница, которая последовательно нарезает блоки xi = (xi,..., xi+m-1)T,

yi = xi+m, i = 1,..., N - m . В результате получаем выборку {Xi, yi}, i = 1,..., N - m , где наблюдение yi играет роль зависимой переменной, а вектор xi выступает в виде вектора независимых переменных. Таким образом, задача прогноза сводится к задаче многомерной регрессии при зависимых наблюдениях. Нетрудно заметить, что в «гусеничных алгоритмах» соседние блоки частично перекрываются, а это приводит к сильной зависимости между ними. Отметим ряд моментов, полученных при анализе условных функционалов для процессов, удовлетворяющих условию слабой зависимости [5]. Условия СП и РСП влекут за собой эргодичность процесса. Свойства непараметрических оценок условных функционалов в случае независимых и слабозависимых наблюдений совпадают - «принцип локальности» («the whitening by windowing principle» Hart (1996)) [5]. «Принцип локальности» имеет место и для случая непараметрических процедур обработки с частично перекрывающимися блоками. В соседних блоках зависимые переменные yi и независимые xi меняются местами. В связи с этим появление даже одного выброса во временном ряду приводит к выбросу как по зависимой переменной, так и к (m-1) выбросам в независимых переменных xi. Следовательно, необходимы робастные алгоритмы регрессии, устойчивые как по оси Х, так и по оси Y. Такие непараметрические алгоритмы имеются [8]. Особую роль для задачи прогноза играет вектор xN-m, который единственен, и только на его основе строится

прогноз. Появление выброса в векторе xN-m приводит к катастрофическим результатам для прогноза. Действительно, весовая функция W(xi / xN-m) играет роль фильтра, вырезающего из временного ряда x1,..., xN-m однородные с xN-m

блоки xi (кластеры). Задача прогноза в этом случае сводится к задаче оценки параметра положения фильтрованных данных yi . Параметр m играет роль размера (интервала зависимости) кластера, а параметр hN определяет ширину пропускания фильтра. Область применения таких алгоритмов связана с некоторыми особенностями исходного процесса X (t), в частности со стационарностью по кластерам. Причем исходный процесс X (t) может состоять из совокупности кластеров с разной волатильностью. Отметим, что с «кластерностью» при изучении случайных процессов столкнулись и при параметрическом подходе при изучении моделей случайных процессов с волатильностью (ARCH, GARCH,...) [1].

2. Устойчивые непараметрические алгоритмы прогноза

В работе [12] предложены методы получения робастных оценок на основе критерия устойчивости. Оказалось, что класс устойчивых оценок может быть синтезирован на основе взвешенного метода максимального правдоподобия (ВММП) [8, 9, 14].

Обозначим через Zi = xi+m - XN+k = yi - XN+k - невязки прогноза и через p(z, 9 / XN-m), P(z, 9 / XN-m) - плотность и функцию распределения невязок в зависимости от реализации вектора XN-m. В данном случае нас интересует неизвестный параметр распределения 9 = (XN+k, S)T . Прогноз XN+k выступает как

параметр сдвига, а оценка параметра масштаба S позволяет оценить вариацию прогноза, то есть

p (x, 9 / XN-m ) =1 h(u), (8)

S

где h(u) - стандартная плотность распределения семейства P(x) и

u = [ x - XN+k ] • S .

Условную М-оценку 9 неизвестного параметра 9 можно определить на основе решения системы эмпирических уравнений вида

jcp,. (t, 9/XN-m Wn (t / XN-m ) = 0, i = 1,2, (9)

где фг- (t, 9 / XN-m) - оценочная функция.

Определим оценочную функцию p(t, 9 / XN-m) в следующем виде [8, 9, 14]:

~ д

ф(^ 9 / Xn-m ) =

g9ln p(^ 9 / Xn-m ) + Р

pl (x, 9 / Xn-m), (10)

где p - параметр, который определяется из условия несмещенности оценки

M[ф(x, 9 / XN-m )] = 0 (11)

и l - параметр, который назовем параметром радикальности оценки.

Выражение (11) определяет ВММП с весами pl (x, 9 / XN-m) : при l = 0 получаем оценки максимального правдоподобия (ОМП); при l = 0,5 - радикальные оценки; при l = 1 - оценки максимальной устойчивости (ОМУ) [12]. Физически роль параметра l сводится к определению степени «мягкого» усечения алгоритма,

настраивая его на вид априорного распределения. Следовательно, варьируя параметром I, можно получать эффективные оценки при локальных отклонениях распределения Р( х, 9 / Хы-т) в классе устойчивых оценок [12].

Для параметров сдвига и масштаба в соответствии с (10) и (12) получаем оценочные функции в виде [9, 14]:

Ф:(х, 9 / XN-т) = к’и (и)к1 -1 (и),

Ф2( х, 9/ XN _т) =

uh’u (u)- 1 1

р1 (х, 9/ / —). (12)

р(х, 9/ Xn _т) (1 +1) _

Система уравнений (9) и оценочные функции (12) определяют алгоритм получения устойчивых оценок параметров сдвига и масштаба в зависимости от параметра радикальности 1.

Рассмотрим в качестве примеров два случая.

■* 1 f (t—X )21

1. p(t, 9 / XN—т) = — exp \-----------N2+k— > - нормальное распределение невя-

v2nS [ 2S J

зок. Из (9), (12) получаем следующую систему уравнений для оценки прогноза.

N—т—к /1 \

Е и eXP \~uf W (Xi / XN—m ) =0; (13)

,■=1 *2

N —т—к / 1 j /• i -j

Е lvU,2 — J+1) eXP (:7U,2 }W (X, / XN—m ) = 0, ()

где W(X, / XN —т ) - определено в (2.10) и и, = [х,+т+ к—1 — XN+к ] ' S^ .

- 1 f It—Xn+J 1

2. p(t, 9 /Xn—т) =— exp ---------------L> - распределение Лапласа. Из (9), (12)

2S I S I

получаем систему уравнений для оценки прогноза:

N—т—к

Е sign(u, )exp{1|u, \}W(X, / Xn—т) = 0; (15)

i=1

N—т—к/ 1 \

Е [l uil — J+1 j eXP {1lui \}W(Xi / XN—т ) = 0. (16)

Можно получить аналогичные уравнения для оценки прогноза и для других

распределений [9]. Сделаем некоторые выводы.

Уравнения (13), (15) определяют:

1) при 1 = 0 - оценки условного максимального правдоподобия (ОМП), которые являются классическими непараметрическими оценками прогноза регрессионного и медианного типов [2, 3, 5, 6, 7];

2) при 1 = 1 - оценки условной максимальной устойчивости прогноза (ОМУ [12]) модального типа [2, 7];

3) при 1 = 0,5 - условные радикальные оценки прогноза [12].

Оценки прогноза ОМП типа (13), (15) не являются робастными по критерию устойчивости (в том числе и оценка медианного типа, которая не справляется с ассиметричными выбросами) [12]. Оценки прогноза модального типа (13), (15) (ОМУ) являются робастными и справляются как с симметричными, так и с асим-

метричными выбросами, но являются низкоэффективными [12]. В связи с этим в [12] и был предложен промежуточный вариант радикальных оценок, у которых эффективность и устойчивость совпадают (они относятся к классу МЭ-оценок на расстоянии Хеллингера). В ВММП радикальные оценки получаются при /=0,5. Алгоритмы вида (13) - (16) были получены в предположении, что вид распределения Р (х, 9 / Хм_т) известен. Следовательно, их можно отнести к частично непараметрическим алгоритмам прогноза. Они являются непараметрическими только по типу априорной информации о виде функции т(Х г _т-1) в модели (1), но не по типу априорной информации о виде функции распределения Е(х, /) е О.

Рассмотрим случай, когда вид распределения Р(х, 9 / Хм_т) неизвестен, и задача относится к классу непараметрических задач.

Пусть Р(х, 9 / Хм_т) относится к непараметрическому классу абсолютно непрерывных унимодальных распределений, симметричных относительно Хм+к. В

этом случае оценка Розенблатта - Парзена для плотности р( х, 9 / Хм _т) запишется в виде

Рм (', 9 / хм _т) = 1|К (■

ХМ+т _ 0,5(^ + 2) £

2 / Хм_т ).

(17)

В соответствии с ВММП для нахождения оценочных функций подставим выражение (17) в (10). В результате получаем оценочные функции

Ф1 (', 9 / Хм_т ) = | К'и (и)^м (2 / Хм_т ) • [Рм V, 9 / Хм_т )] _' ,

Ф2 (', 9 / ^^м_т ) = {/ иКи (и)^м (2 / ;Тм_т ) • [Рм а, 9/ ^Тм_т )]“' _ ^

Xр1м (Г, 9/ Xм_т ), (18)

где и = [Xм+т _ 0,5(/ + 2)] • £_ .

Подставляя (18) в (9), получаем систему оценочных уравнений для нахождения непараметрических оценок прогноза Хм+к и параметра масштаба £ :

|фг (/, 9/Хм_п )^ а/Хм_п) = 0, г = 1,2. (19)

Пусть

м+т

= 0,5 •[ хг-

Оценочные уравнения (19) примут следующий вид:

м_т_к м_т_к

X X К' (и„ )Ж(ху / хм_т )

г=1

у=1

м _т_к

X К(и„ )Ж(ху / хм_т )

У=1

I _1

м—т_к м_т_к

X игу • К' (игу )Ж(^ / хм_т ) '

V=1

м _т _к

X К (Uгv )Ж (хv / хм

|_ v=1

1

7+1

м _т_ к

X К(ип )Ж(Xv / хм_т )

Ж(х / хм_т ) = 0.

(20)

Оценочные уравнения (20) выглядят достаточно сложно, особенно если подставить выражения для ядерных функций К(иіу) и Ж(хі / хы-т). Однако при 1=0 (ОМП) и применении «идеального окна» К (и) получаем хорошо известные в непараметрической статистике алгоритмы оценки сдвига и масштаба на полусуммах Уолша.

Рассмотрим влияние одиночных выбросов на работу непараметрического алгоритма (20). Пусть выбросом является одно из наблюдений х1,...,xN-m . Если выброс входит в блок хі, то фильтр Ж(хі / хы-т) приписывает этому блоку малый вес или просто не пропускает данный блок. Непараметрический алгоритм ведет себя устойчиво. Если выброс пришелся на значение уі, то алгоритмы типа «радикального» или ОМУ решают задачу робастного оценивания прогноза, в том числе и для ассиметричных распределений выбросов [8, 12]. Если выбросы отсутствуют, то при I = 0 алгоритм (20) сходится к оптимальному алгоритму. Таким образом, ВММП позволяет, применяя адаптивную настройку параметра радикальности 0 < I < 1, получать эффективные непараметрические оценки прогноза в неопределенной ситуации с выбросами [8].

Пусть выброс присутствует в векторе хы-т. В этом случае фильтр

Ж(хі /хы-т) не найдет во временном ряду х1,...,хы-т ни одного вектора хі, однородного с вектором хы-т. На выходе непараметрического алгоритма прогноза получим Хы+к = 0. Данная информация важна, например, в задачах определения момента перескока процесса на новый уровень. Она сигнализирует о том, что в векторе хы т присутствует выброс. Обычные методы борьбы с выбросами, связанные с исключением данного вектора или занижением его роли с помощью специальных весовых функций, в этом случае не пригодны. В дальнейшем рассмотрим два алгоритма, позволяющие учесть данную ситуацию.

Теоретическое исследование алгоритмов вида (20) выходит за рамки данной работы и связано с исследованием непараметрических оценок параметров, заданных в неявном виде через обобщенные условные функционалы (обобщенные условные М-оценки)

|ф(Е(х, 0 /ґ)Е’(х, 0 /ґ),...)йЕ(х, 0 /ґ = 0 .

3. Адаптивные алгоритмы

Эффективность работы алгоритмов вида (20) определяется набором параметров (т, кы) для настройки фильтра Ж(хі / хы-т) и параметром радикальности I для настройки робастности алгоритма. Для оптимизации по этим трем параметрам возможны два подхода. Первый основан на теоретических исследованиях вариации оценки прогноза и ее оптимизации по данным параметрам. Второй подход, который в основном и применяется, основан на минимизации наблюдаемой суммы квадратов невязок Q(m, км, I) = (N - т)-1 ^ (хі+т - Хы+к )2 в режиме «скользящего экзамена».

Уже предварительные результаты показывают, что задача минимизации Q(m, ^, I) хорошо сегментируется, хотя все три параметра связаны. В режиме

«скользящего экзамена» сначала настраиваются параметры (т, Ны) фильтра Ш(хі / хм_т), которые слабо зависят от выбросов внутри временного ряда х1,..., хы_т , а затем параметр радикальности I.

Как было отмечено в пункте 2.2, непараметрические алгоритмы прогноза не спасают от выбросов в векторе хы_т. Для отбраковки выбросов в хы_т можно предложить два метода: прореживания и прореживания с восстановлением. Предположим, к=1 и выброс приходится на значение хы .

Прореживание. Исключаем наблюдение хм из вектора хм_т и будем проводить прогноз по вектору хм_т_1 = (хм_1,...,хм_т_1)Т на два шага вперед. Если выброс находится на у-м месте вектора хм_т, то, последовательно удаляя и возвращая в хм_т по одному наблюдению, обнаружим выброс и найдем оценку прогноза. Во временных рядах, в связи с сильной зависимостью между соседними наблюдениями, могут возникать пачки выбросов. В этом случае исключаем пары, тройки, ... соседних наблюдений в векторе хы_т.

Прореживание с восстановлением. Сначала с помощью прореживания выявляется выброс, например хы . Затем наблюдение хы с помощью непараметрического алгоритма прогноза восстанавливается хы и прогноз ведется по восстановленному вектору хы _т = (% , % ^..^ % _т )Т .

Проведены исследования алгоритмов (13), (14) и (19) с нормальными ядрами в предположении, что выбросов в %-т нет. В первую очередь интересовало поведение Q(m, кы, I) в зависимости от параметра т (размер кластера) и тесно связанного с ним параметра кы (ширина фильтра). Затем, как и в работе [8], производилась адаптация по параметру радикальности I. Модель процесса Х1У) = 5ш(ю- t + ф) + 0,15е, е - белый шум, ф - равномерна на [0,2п] и модель биений X2 У) = • t + ф) БШ(ю2 • I) + е .

На рис. 1 - 3 представлена зависимость Q(m,кы,I) от т для X1(t): рис. 1 -без шума; рис. 2 - присутствует шум е; рис. 3 - присутствует шум е и единичный выброс.

4. Моделирование

Є-10-2

Є-10-2

Є

3.5 3

2.5 2

1.5

6

7

5

0,23

4

1

0,5

2 4 6 8 10 т

3

2 4 6 8 10 т

2 4 6 8 10 т

Рис. 1

Рис. 2

Рис. 3

Как видно из рис. 1, уже при m = 2 имеет место сингулярный прогноз. Наличие белого шума (рис. 2) приводит к регулярному прогнозу и m увеличивается до 8. Наличие белого шума и выброса (рис. 3) приводит к регулярному прогнозу, причем m не увеличивается, а адаптация по l приводит к ликвидации влияния выброса. Зависимость Q(m, hN, l) от l качественно имеет такой же вид, как и в [8], поэтому не приводится. Аналогичная картина наблюдалась и для процесса биений X 2(0.

Заключение

В работе на основе ВММП синтезирован ряд новых робастных непараметрических алгоритмов прогноза (алгоритмы (10), (13), (15), (20)). Данные алгоритмы включают как частный случай непараметрические алгоритмы [2, 3, 5 - 7].

Показано, что введенные ранее непараметрические алгоритмы прогноза [2, 3, 5 - 7] не являются полностью робастными и непараметрическими.

Предложены адаптивные методы настройки по параметрам m, hN, l, которые оптимизируют работу алгоритмов как по исходному распределению, так и по распределению выбросов.

Проведено моделирование, которое показало работоспособность и эффективность алгоритмов.

В теоретическом плане алгоритмы вида (20) ставят вопрос исследования непараметрических оценок параметров, заданных в неявном виде через обобщенные условные функционалы (обобщенные условные М-оценки).

ЛИТЕРАТУРА

1. Ширяев А.Н. Основы стохастической финансовой математики. T. 1. М.: Фазис, 1998. 488 с.

2. Siegfried Heiler. A Survey on Nonparametric Time Series Analysis. 1999. 49 p. URL: http:// www.ub.uni-konstanz.de/kops/volltexte/1999/316/

3. Neumann M.H., Kreiss J.P. Regression-Type inference in nonparametric autoregression // Ann. Statis. 1998. V. 26. №. 4. P. 1570-1613.

4. Мартин Р.Д. Устойчивый авторегрессионный анализ временных рядов // Устойчивые статистические методы оценки данных. М.: Машиностроение. C. 121-146.

5. Васильев В.А., Добровидов А.В., Кошкин Г.М. Непараметрическое оценивание функционалов от распределений стационарных последовательностей. М.: Наука, 2004. 510 с.

6. Симахин В.А. Непараметрическое прогнозирование случайных процессов // Тез. докл. I областной научно-практической конференции по надежности научно-технических прогнозов. Новосибирск, 1978.

7. Рымар Т.Н., Симахин В.А. Непараметрическое прогнозирование стационарных случайных процессов // Тез. докл. зональной научно-технической конференции «Датчики и средства первичной обработки информации». Курган, 1990. C. 102-104.

8. Simahin V.A. Nonparametric robust regression estimate // Proc. SPIE. 2006. V. 6522. P. 130-139.

9. Симахин В.А. Непараметрическая статистика. Ч. II. Теория оценок. Курган: Изд-во КГУ, 2004. 163 с.

10. Кошкин Г.М. Об одном подходе к оцениванию переходной функции распределения и моментов для некоторых марковских процессов // Математическая статистика и ее приложения. Томск, 1976. C. 53-65 .

11. Roussas G.G. Nonparametric estimation of the transition distribution function of a Markov process // Ann. Math. Statist. 1969. V. 40. №>. 4. P. 1386-1400.

12. Шурыгин А.М. Прикладная статистика. Робастность. Оценивание. Прогноз. М.: Финансы и статистика, 2000. 223 с.

13. Симахин В.А. Непараметрическая статистика. Ч. I. Теория оценок. Курган: Изд-во КГУ, 2004. 216 с.

14. Симахин В.А. Взвешенный метод максимального правдоподобия II Материалы IX Международной науч.-технич. конф. «Кибернетика и высокие технологии XXI века». Воронеж, 200В. T. 2. C. 661-672.

Симахин Валерий Ананьевич

Курганский государственный университет (г. Курган).

E-mail: [email protected] Поступила в редакцию 16 сентября 2010 г.