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

УДК 519.242.33

Ю. Н. Орлов, А. Д. Босов

Московский физико-технический институт (государственный университет)

Кинетико-гидродинамический подход к прогнозированию нестационарных временных рядов на основе уравнения Фоккера^Планка

Выводится цепочка кинетических уравнений эволюции многомерных выборочных функций распределения нестационарного временного ряда. Строится система уравнений гидродинамики, замыкание которой на каком-либо моменте распределения приводит к модели временного ряда в виде порождающей его нелинейной динамической системы.

Ключевые слова: выборочная плотность функции распределения, уравнение Лиувилля, уравнение Фоккера—Планка, нестационарный временной ряд.

1. Введение

При анализе нестационарных временных рядов одной из центральных задач является оптимизация объема выборки для вычисления текущих статистик ряда. Необходимость оптимизации связана с тем, что малый объем статистически мало репрезентативен, а большой объем данных не дает нужной точности в силу того, что текущее распределение изменилось. Для практического решения этой задачи в работах [1,2] была введена индикативная статистка, названная авторами «горизонтным рядом» и показывающая, каким объемом лучше всего, т.е. с минимальной ошибкой в среднем квадратичном, определяется текущая выборочная плотность функции распределения (далее ВПФР).

Пусть х(Ь) есть значение случайной величины £ в момент времени £, и /т(х,£) есть ВПФР, построенная по выборке объема Т в момент времени Ь. Говоря о случайной величине £, мы предполагаем, что, по крайней мере, теоретически в каждый момент времени существует некоторая генеральная совокупность, определяющая вероятность появления того или иного значения ж. Для простоты ползаем, что время £ дискретно, ряд х(Ь) эквидистантный с единичным шагом, причем все значения равномерно ограничены по времени. В этих условиях без ограничения общности можно считать, что V Ь х(Ь) € [0; 1]. Основное утверждение, показывающее возможность равномерного ограничения по времени для оценки неточности в определении ВПФР, состоит в том, что для функционала V(Т, т; Ь) близости двух ВПФР

|/т(х,г + т) - /т(М)| Лх (1)

о

V (Т,т; 1) = имеет место оценка [1]

0 < V(Т,т; г) < шт (2т/Т; 2). (2)

[0; 1]

ВПФР в граничных точках полагаются равными нулю. Пределы интегрирования далее для краткости не указываются.

Определение 1. Горизонтным рядом к(Ь) для ряда х(Ь) при сдвиге на промежуток т называется такой минимальный объем выборки к (Ь, т; в), что при всех Т ^ Н(Ъ,т; е) выполнено условие V(Т, т; Ь) ^ е.

Таким образом, горизонтный ряд - это двупараметрическая статистика, которая характеризует взаимозависимость двух выборок. Оптимальный объем выборки находится из распределения горизонтного ряда, если последнее квазистационарно на горизонте прогнозирования т, т.е. если различия в распределении горизонтной статистики составляют величину о(е). Примеры нестационарных рядов с таким свойством приведены в [3].

Оптимизация объема выборки приводит к тому, что распределение /т(х, Ь) на горизонте т в первом приближении (с точностью е) можно считать стационарным. В следующем приближении можно построить уравнение, описывающее эволюцию ВПФР на этом горизонте. В силу нормировки ВПФР на единицу это уравнение должно иметь вид закона сохранения, т.е. быть аналогом уравнения Лиувилля в классической статистической механике. В [3] показано, что можно придать корректный смысл эмпирическому уравнению Лиувилля:

где локальная скорость ит(х,Ь), называемая в [3] эмпирической лиувиллевой скоростью (ЭЛС), определяется через двумерную ВПФР Р (х,х,Ь) значений ряда и его приращений (далее указание на объем выборки Т для краткости опускаем) по формуле

Производные по времени в (3) и (4) трактуются как разностные производные с шагом 1,

[0; 1]

Эта трактовка оправданна, если Т >> 1. Если удастся каким-либо образом определить ЭЛС, то из (3) можно будет найти и ВПФР. Прогноз самого ряда строится как среднее значение по найденной ВПФР.

В настоящей работе предложен метод определения ЭЛС с помощью системы уравнений для эволюции моментов ВПФР /(х,Ь), которые замыкаются на третьем уравнении. В гидродинамической аналогии это уравнение энергии. Замыкание происходит в духе «наивного прогноза», когда момент третьего порядка считается постоянным. Корректность такой гипотезы проверяется постфактум, т.е. модель признается приемлемой для данного временного ряда, если изменение «постоянной» величины на горизонте т можно считать величиной порядка о(е). На практике обычно полагают, что е = 0, 05, а о(е) ^ 0, 01.

Прогнозная модель на основе уравнения Лиувилля (3) была предложена в [3]. Здесь будет построена модель, основанная на эволюционном уравнении Фоккера—Планка, которое моделирует не только снос плотности вероятности со скоростью и(х,£), но и эффект случайного блуждания, оцениваемого по эмпирическим данным с предыдущего шага по времени. Эта модель оказалась более точной в плане прогнозирования временного ряда, чем модель, основанная на уравнении Лиувилля.

2. Цепочка кинетических уравнений для ВПФР

Опишем идеологию вывода кинетических уравнений для эволюции ВПФР, следуя подходу А.А. Власова, изложенному им в [4]. Предположим, что в нашем распоряжении имеется вся совокупность величин х(Ь), так что можно построить совместные «генеральные» плотности функции распределения (ПФР) /п (х, х, ..., х(п); любой размерности. Мера

представляет число кинематических состояний частиц в границах указанных промежутков. Согласно формулам теории вероятностей для плотностей совместных многомерных распределений, имеет место зацепление ПФР в виде

д/т (х,г) + д_

дЪ дх

+ — ит (х,г)/т (х,£) = 0,

(3)

и(х,£)/(х,£) = хР(х,х,£)йх, /(х,£) = Р (х,х,1) .

(4)

(5)

/п{.

X X X

{п); *) = /п+1 (:

X X X

(п), х(п+1),^1 с1х(п+1)

(6)

Нетривиальным является вопрос о том, что представляет собой математический объект, получающийся из (5) предельным переходом п —^ В [4] этот вопрос не обсуждался;

видимо, предполагалось, что смысл объекта

ГО

dГ = Р (х, х, х, ..., х(к\ ...; ^ ^ dx(k') (7)

к=0

интуитивно понятен. Однако далее для плотности меры (7) будет записано некое эволюционное уравнение, поэтому требуется определить ее более корректно. Такое определение дано в недавней работе В.В. Козлова и О.Г. Смолянова [5], где построена также и общая схема перехода от уравнения Лиувилля относительно мер на бесконечномерных пространствах к равносильной этому уравнению бесконечной системе уравнений относительно функций распределения на конечномерных пространствах. Мы не будем углубляться в теорию, развитую в [5], поскольку наша задача формулируется в терминах конечномерных пространств, но упомянуть об этой проблеме было необходимо для целостности концепции. Таким образом, можно считать, что уравнение Лиувилля для Р, выражающее закон сохранения нормировки, корректно, как и представления проекций на конечномерные подпространства:

I + Е?=оА**» №М)Р) =0; т

( о )

/п (х, X, ..., х(п); *) =/ Р(х,х ,х, ...,х(к\ ...; {) Пь=п+1 Лх(к).

После этого, считая, что при стремлении скоростей и ускорений всех порядков к бесконечности ПФР стремится к нулю, из (8) и из соотношений (6) получаем уравнения эволюции для п-частичных ПФР:

^ + А1уж § х/1 (х, х; ^ dx = 0 ;

дь{х,х,г) + а\ухх/1 (х, х; ^ + ^у* §х/2 (х, х, х; $ dx = 0 ; ... ;

ди(х,х,^,х( ),£) + ^_о х(к+1)^п (х, ±, ..., х(п); ^ +

+А!уж(п) §х(а+1')/п+1 (х, х, ..., х(п+1'); £) dx(n+1') = 0 ; ... .

(9)

Таким образом, уравнение эволюции для каждой частичной ПФР не является замкнутым. Оно связывает между собой концентрации и потоки, причем последние выражаются через ПФР от большего числа переменных. Если на каком-либо этапе оборвать зацепление, т.е. задать последний интеграл в (9) независимо, то в пределе нулевых центральных моментов всех порядков относительно координат, скоростей и ускорений всех порядков будет получаться та или иная динамическая система [3].

Если замкнуть цепочку (9) на первом уравнении, считая интеграл в нем равным нулю, получится уже упомянутое выше в п. 1 стационарное решение. В следующем приближении интеграл по учитывается, а обрыв происходит на втором уравнении, где интеграл по dx считается равным нулю и т.д. Так, при замыкании на втором уравнении система уравнений для /о(х,Ь) и $1(х,х,£) имеет вид

+ I иШШх,0 = 0, +х =0; (10)

и(х, 1)/о(х, ^ = $ х/1(х, х, t)dx , /о(х, ^ = / /1 (х, х, I) dx .

Заметим теперь, что в системе (10) скорость и(х,£) может быть любой, лишь бы при решении первого уравнения не появились бы отрицательные вероятности. Статистический смысл ЭЛС состоит в том, что ее среднее значение равно производной среднего выборочного значения ряда [2]:

d

dt

х/о(х, t)dx =

д/о(х,г)

х------------dx = —

о1

д

х—и(х,^)/0(х,1)йх =

и(х,1)/0 (x,t)dx = и (I).

(11)

Действительно, эта же величина равна изменению эмпирического первого момента за один шаг по времени:

1 1 1 т{г+1)-т{г) = - ^2 х{Ъ—Т+к+1)--^2 %(І-Т+к) = - ^2 %(І-Т+к)

к=1 к=1 к=1

х/\(х, X, і)д,х,

а последнее выражение равно (11) в силу определения ЭЛС (4).

Очевидно, наблюдаемое изменение выборочного среднего значения временного ряда может быть получено при различных ЭЛС и(х,£), которые в этом смысле определены С ТОЧНОСТЬЮ до выражения вида \(Ь)д/о(х,1)/дх, где Х(Ь) — произвольная функция времени. Это наблюдение позволяет расширить класс кинетических моделей для описания эволюции ВПФР.

3. Эмпирическое уравнение Фоккера—Планка

Цепочка равенств (11) показывает, что эволюция первого момента ВПФР в силу уравнения Лиувилля (3) совпадает с производной по времени среднего выборочного значения, вычисленного непосредственно по временному ряду. Однако для высших моментов такой инвариантности нет. Тем не менее можно так «подправить» ЭЛС без нарушения нормировки ВПФР, чтобы первые два эмпирических момента ВПФР изменялись бы в соответствии с кинетическим уравнением эволюции самой ВПФР.

Теорема 1. Пусть т(Ь) и д(Ь) — первый и второй выборочные моменты нестационарного временного ряда, /0(х,Ь) — его ВПФР, и(х,Ь) — ЭЛС, а и(Ь) — среднее значение

ЭЛС. Тогда, эволюция величин т(Ь) и д{Ъ) в силу уравнения Фоккера Планка:

д/0 д(и/0) ,,,д2/0

■Ж + ^ — а? = °- <12>

где

\(Ъ) = 1 д(1) — соу(х, и) — т(1)и(I), (13)

совпадает с их изменениями, вычисленными непосредственно по выборке.

§

т(Ь) = |х/0(х,1)(1х. Рассмотрим далее производную (1д/(И в силу уравнения (12):

2 д/о,

х —— ах = оі

(-

2 1 д(и!о) , ^92fо\ ^

х\-------------------------------я-+ Лх-

ох ох2 I

После двукратного интегрирования по частям с учетом того, что на границах промежутка ВПФР равна нулю, получаем

{и/о - Х1Б)

и!о — ) (ІХ = 2соу(х, и) + 2т,и + 2\,

где со\(х, и) обозначает ковариацию величин х и и(х,ї). Полученное выражение совпадает с выражением (13). Теорема 1 доказана.

Таким образом, величина \(Ь), представляющая собой коэффициент диффузии, равна половинной невязке между производной по времени второго выборочного момента ВПФР и его эволюции в силу уравнения Лиувилля (3). Уравнение (13) будем называть эмпирическим уравнением Фоккера—Планка.

Может показаться, что замена уравнения Лиувилля (3) на уравнение Фокера—Планка (13) ничем не лучше, поскольку в момент времени і, производная д(і) не известна, так как для ее независимого определения нужно знать значение ряда в момент і + 1. Тем не менее модель (13) можно использовать на практике после небольшой коррекции, состоящей в

х

замене А(£) на Х(Ь—1). На следующем шаге по времени этот коэффициент, естественно, пересчитывается с учетом нового значения х(Ь+1). Ошибка, которая вносится при такой коррекции, часто оказывается величиной порядка о(е). Это связано с тем, что для е-стационарных ВПФР, т.е. таких, для которых выполнено условие V(Т, 1;Ь) = ||/(х,Ь) — /(х,Ь + 1)|| ^ £, изменение второго момента на [0; 1] сравнительно мало. Действительно, пусть N есть число промежутков разбиения отрезка [0; 1] при практическом вычислении интеграла нормы V(Т,т; Ь), и /(г, Ь) есть значение ВПФР в г-м промежутке. Тогда из того, что

N

V г \/(г^) — + 1)| < ^\1(г^) — /(г,г + 1)| < е,

г=1

следует оценка

\д(Ь + 1) — д(Ь)\ = |Г х2/о(х, t + 1)с!х — Гх2/0(х, 1)с!х\ ^

1 (14) х<2 \/о(х,Ъ + 1) — !о(х,^) \ Лх ^ £/0 х2(1х = § .

В то же время, как нетрудно убедиться, невязка производной второго момента ВПФР с вычисленной производной по уравнению Лиувилля составляет величину порядка е. Таким образом, даже скорректированное уравнение Фоккера—Планка дает более точное приближение к эволюции моментов ВПФР, чем уравнение Лиувилля.

4. Система уравнений гидродинамического типа и ее замыкание

Кинетическая модель эволюции ВПФР /о(х,Ь), основанная на уравнении Фоккера— Планка (12) и замкнутая на уравнении для двумерной ВПФР /о(х, х, £) так, как это сделано в (10), имеет вид

^ =0, = 0; (15)

и(х, {)/о(х, ^ = $ х/0(х, х, 1)д,:Ъ , /0(х, ^ = / /0 (х, х, {) .

Заметим, что сама ВПФР /о(х,Ь) построена по соответствующему оптимальному объему выборки при прогнозировании на 1 шаг вперед, как описано в п. 1. Вычисление же совместной ВПФР /о(х,х,Ь) требует гораздо большего объема машинной памяти, не говоря уже о работе со следующим приближением (обрыве на третьем уравнении цепочки (2.5) и т.д.). Поэтому практически важной задачей является создание эволюционной модели непосредственно в терминах одномерной ВПФР. Из (15) следует, что для этого надо определить ЭЛС и(х, Ь). Это можно сделать, либо найдя /0(х,х,Ь) из оборванного уравнения Лиувилля для парной ВПФР, что требует больших вычислительных затрат, либо получить уравнение для самой и(х,Ь).

Используя второе и третье уравнения системы (15), получаем

9 (и/о)

дИ

. ЩхЛ л. = —

о1

^2 д/о(х,х,г) л±= д(е(х,г)/о(х,г)) дх дх ’

где е(х,£) представляет, как это видно из (16), плотность второго момента по скорости для ВПФР /о(х,х,Ь). Таким образом, эволюция и(х, Ь) (т.е. первого момента /0(х,х,Ь)) определяется вторым моментом. Возникает эволюционная цепочка для моментов, решение которой эквивалентно построению эволюции самой ВПФР в предположении, что проблема моментов однозначно разрешима. Но как следует из оценки (14), при нахождении приближенной эволюции ВПФР нет необходимости строить уравнения эволюции для моментов всех порядков, поскольку, чем выше порядок момента, тем меньше его вклад в эволюцию ВПФР. Следовательно, моментную цепочку также можно оборвать на каком-либо порядке,

считая интеграл в правой части известным. Тогда получается замкнутая система уравнений гидродинамического типа. Например, при замыкании на третьем уравнении (т.е. на эволюции момента е(х,£)) получаем

+ Ши(х,1)/о(х,1) — Щ)= 0 , (17)

Ли(х, Ь)/о(х, Ь) = — ^е(х, Ь)/о(х, г), ^е(х, Ь)/о(х, г) = — Ц:д(х, Ь)/о(х, Ь),

где д(х,£) = /х3/о(х,х,£)йх считается известной функцией. Например, эту функцию, как и величину \(Ь), можно считать известной с предыдущего шага по времени.

5. Пример модели прогнозирования нестационарного временного ряда

Опишем кратко результаты применения модели (17) к задаче прогнозирования конкретного временного ряда, образованного суточными данными с(Ь) курса ЦБ рубль/доллар за период с 2005 по 2010 гг. по данным [6]. Обозначим

М = шах (с(1 + 1) — с(Ъ)) , т = шт (с(1 + 1) — с(1)) (18)

и рассмотрим нормированные приросты курса:

т ф —. (и)

Ряд х(Ь) существенно нестационарный, что может быть проверено стандартными методами, например по критерию Колмогорова—Смирнова [7]. Так если две выборки объемом в 500 данных, расположенные встык, взяты из одной генеральной совокупности, то они отличаются в норме непрерывных функций не более чем на 0,1 с вероятностью 0,9. Однако оказалось, что из 600 таких выборок, построенных скользящим образом по рассматриваемым данным, лишь 28 (менее 0,05) отличаются менее чем на 0,1, тогда как в среднем

отличие между ними составило 0,23, а 0,9-квантиль распределения расстояний между выборками лежит от 0,11 до 0,35.

Исходный ряд с(Ь) является настолько нерегулярно меняющимся, что наилучшим его стандартным прогнозом (т.е. адаптивным, регрессионным или авторегрессионным, скользящим средним и т.п.) оказывается наивный (завтра будет как сегодня), ошибка которого равна среднеквадратичному относительному приросту и составляет для изучаемого фрагмента примерно 0,008. Ошибка вышеперечисленных стандартных прогнозов по оптимальному для них объему данных заключена в диапазоне от 0,009 до 0,011, т.е. все они работают хуже наивного прогноза.

Рекуррентно разрешимая модель для построения ВПФР /о(х,Ь), имеющая вид нелинейной динамической системы, получается из (17) путем дискретизации области значений х(Ь). Величина ячейки при этом должна иметь порядок о(е). Например, при е = 0, 05 можно взять 100 ячеек равномерного разбиения. Тогда, обозначая тильдой рассчитываемые величины, получаем в праворазностном координатном шаблоне следующую систему:

{ : / (х,1 + 1) = /(х, I) + /(х, {)й(х, I) — /(х + 1,{)й(х + 1,1) +

+(/ (г + 1,1) — 2/(г, I) + / (г — 1,1)) ;

и : (й(х, I) — и(х, £ — 1)) /(х, £ — 1) =

= и(х, £ — 1) (и(х + 1,1 — 1)/(х + 1,1 — 1) — и(х, £ — 1)/(х, £ — 1)) +

+ ё(х, {)/(х, £ — 1) — ё(х + 1,1)/(х + 1,1 — 1) ; е : (ё(х, I) — е(х, £ — 1)) /(х, £ — 1) =

= е(х, £ — 1) (и(х + 1,1 — 1)/(х + 1,1 — 1) — и(х, £ — 1)/(х, £ — 1)) +

+ д(х, £ — 1)/(х, £ — 1) — д(х + 1,1 — 1)/(х + 1,1 — 1) .

Анализ горизонтной статистики для ряда х(Ь) величин, определенных в (19), дает оптимальный объем для прогнозирования на 1 шаг вперед, равный в среднем Торъ = 290 данных. Эволюция такой ВПФР по системе уравнений (20) отличается от фактической на

0.003. что более чем в 2 раза меньше, чем верхняя грань ошибки, равная, согласно (2), 2/Торх. Прогноз самого ряда как среднего значения прогнозной ВПФР имеет относительную ошибку отклонения от факта, равную в среднем квадратичном по рассматриваемому фрагменту данных 0,004. Это в 2 раза лучше, чем ошибка наивного прогноза.

Таким образом, существуют примеры рядов, когда предложенная модель имеет заметные преимущества в точности перед стандартными методами прогнозирования, корректно применяемыми лишь для стационарных рядов. Дальнейшее уточнение модели возможно на пути более полного учета высших моментов ВПФР совместных приращений высших порядков для случайных величин. Ограничением здесь является лишь производительность вычислительных комплексов, поскольку с увеличением размерности объем анализируемых данных стремительно возрастает. Например, анализ одной только горизонтной статистики для трехмерной задачи (с учетом разностных ускорений) на том же объеме данных (порядка 1000) требует использования массива 1016 точек. Заметим, что объем выборки в тысячу данных сравнительно мал, на практике более востребованы задачи по анализу выборок объемом несколько миллионов данных. Разработка соответствующих численных алгоритмов представляется актуальной и практически важной задачей.

Литература

1. Орлов Ю.Н., Осминин К.П. Методика определения оптимального объема выборки для прогнозирования нестационарного временного ряда // Информационные технологии и вычислительные системы. — 2008. — № 3. — С. 3-13.

2. Орлов Ю.Н., Осминин К. П. Построение выборочной функции распределения для нестационарного временного ряда // Мат. мод. — 2008. — №9. — С. 23-33.

3. Орлов Ю.Н., Осминин К.П. Нестационарные временные ряды. - М.: Эдиториал УРСС, Книжный дом ЛИБРОКОМ, 2011.

4. Власов А.А. Статистические функции распределения. - М.: Наука, 1966.

5. Козлов В.В., Смоляное О.Г. Бесконечномерные уравнения Лиувилля относительно мер // Доклады РАН. - 2010. - Т. 432. - № 1. - С. 28-32.

6. Данные обменных курсов валют // www.finam.ru

7. Гнеденко Б.В. Курс теории вероятностей. - М.: Физматлит, 1961.

Поступим в редакцию 22.02.2012.