CC BY
88
УДК 519.6+004.94 ББК 22.193
МОДЕЛИРОВАНИЕ РЕШЕНИЙ СТОХАСТИЧЕСКИХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С МАРКОВСКИМИ ПЕРЕКЛЮЧЕНИЯМИ
Черных Н. В.1
(Арзамасский политехнический институт (филиал) Нижегородского государственного технического университета им Р.Е. Алексеева)
Рассматриваются математические модели сложных систем в виде стохастических дифференциальных уравнений с марковскими переключениями диффузионной составляющей, зависящими от фазового состояния. Проводится моделирование решений таких уравнений в среде Scilab с использованием схем Эйлера, Мильштейна, Тейлора, сравнивается качество аппроксимаций.
Ключевые слова: стохастические системы, марковские переключения, схема Эйлера, схема Мильштейна, схема Тейлора для стохастических систем, сходимость, устойчивость, погрешность, доверительный интервал.
1. Введение
Вследствие широкого диапазона применений, гибридные системы, также известные как системы с переключаемой диффузией, стали более популярными в последние годы и привлекают повышенное внимание, особенно в области контроля и оптимизации. Как отмечают в [25] G. Gerorge Yin, Chao Zhu, одной из причин популярности систем с переключаемой диффузией
1 Надежда Валентиновна Черных, аспирантка (nadezdacher@mail. ru).
является то, что много систем реального мира в новую эру требуют построения сложных моделей; кроме того, увеличены требования к моделированию больших систем, проектированию оптимальных средств управления. Объединяя непрерывную динамику и дискретные события, гибридные системы способны описывать сложные системы, их изменчивость и хаотичность в окружающей среде.
Гибридные системы рассматривали для моделирования систем электроэнергии Willsky и Levy, для контроля солнечного теплового центрального приемника - Sworder и Rogers и многие другие. Еще в 1987 г. Athans в работе «Command and control (C2) theory: a challenge to control science» предположил, что гибридные системы станут основной структурой в изложении и решении связанных с контролем проблем в крупных системах управления и коммуникаций [13]. Гибридные системы с марковским режимом переключений использовались для моделирования многих физических систем, которые могут испытывать частые непредсказуемые структурные изменения, вызванные явлениями типа повреждения или отказа компонент и резких волнений окружающей среды. Ссылки на эти и другие работы можно найти в [13, 15 , 25].
Процессы с механизмом переключений используются также для моделирования производственных процессов, описания функционирования биржи, для сбора данных об отдельных изменениях в финансовой сфере и страховании, а также решения стохастических проблем оптимизации в беспроводных каналах связи и сетях.
Современные финансовые рынки являются сложными системами, состояние которых может резко изменяться во времени с восходящими и нисходящими тенденциями. Их изучению и моделированию посвящено много исследований. Например, в [17] авторы исследуют драматические события экономических кризисов, сравнивая мировой финансовый кризис последних лет с депрессией 1930-х. В [6] авторы отмечают, что модель марковской переключаемой диффузии с непрерывным временем широко используется в финансовой математике, и предлагают свой
метод оценки параметров для реального соответствия ситуации на мировом рынке. В последние годы страхование жизни становятся все более популярным в крупных страховых и финансовых рынках во всем мире. Заключая контракт, держатели страхового полиса платят ежегодные взносы страховщику, который управляет капиталом. Одна ключевая особенность этих инвестиционных планов - разделение прибыли от инвестиционного портфеля между держателями страховых полисов и страховщиком. В [21] предлагается модель, позволяющая оценить результаты участия в таких проектах, описанная переключаемой диффузией под управлением марковской цепи. В [6, 7, 11, 17, 21] можно найти ссылки на множество других работ в этой сфере.
Моделирование и имитация биохимических систем - важные задачи, потому что они могут обеспечить понимание сложной системы в случаях, когда традиционное экспериментирование дорого или невозможно. Биологические системы - часто смеси непрерывных и дискретных процессов со сложной динамикой, в которых разъединение и изучение отдельных компонент может пролить свет на функции всей системы. Моделирование таких систем - интересная, перспективная, но сложная задача. Стохастические гибридные системы - идеальное средство для моделирования биохимических систем. В работе [20] авторы развивают метод моделирования стохастической гибридной системы, а также представляют социологическое исследование системы баланса воды/электролита в людях. Там же можно найти ссылки на другие работы.
В [5] авторы проводят вероятностный анализ пространственно-временного распространения гриппа в Германии, также используя диффузионной моделью под управлением дискретной цепи Маркова. Грипп - одна из самых распространенных и серьезных болезней во всем мире, поэтому важно понимать характер распространения болезни. Вагоуап, Rvachev и Иванников (1977) были первыми, кто смоделировал пространственное распространение гриппа, рассматривая сеть транспортировки из определенного региона на примере системы железнодорожного транспорта в СССР.
Один из важных классов гибридных систем - стохастические дифференциальные уравнения с марковскими переключениями (SDEwMSs). Большинство SDEwMSs не имеет явных решений, поэтому важно иметь численные решения. В 2004 г. C. Yuan, X. Mao отмечали в [28], что «нет никаких численных методов для SDEwMSs, хотя методы решения стохастических дифференциальных уравнений (SDEs) хорошо изучены», и предлагали свой метод решения подобных уравнений.
В последующие годы численное решение SDEwMSs изучалось интенсивно разными авторами. Ссылки на такие работы можно найти в [12, 14, 22, 26, 28]. В 2010 г. в свой книге [25] G. Gerorge Yin, Chao Zhu отмечали, что, несмотря на то, что численные методы для SDEwMSs были рассмотрены многими исследователями, меньше известны методы моделирования процессов с переключениями, зависящими от состояния, т. е. изучение такого механизма переключений находится все еще в периоде становления, так как в нем есть свои трудности.
Среди работ последних лет по развитию методов решения SDEwMSs и оценке погрешности можно отметить [4, 9, 10, 16,
18, 19, 23].
Данная работа продолжает исследование, начатое в [3, 25, 27] по моделированию решений SDEwMSs с переключениями, зависящими от фазового состояния.
Многие авторы работ по численному моделированию отмечают, что сходимость численной схемы не является гарантией того, что схема может использоваться эффективно в практическом моделировании. Неустойчивость схемы может сделать ее бесполезной для решения практических задач.
В данной работе будем сравнивать качество аппроксимаций, полученных с использованием известных численных схем Эйлера, Мильштейна, Тейлора, адаптированных для моделирования переключаемых диффузионных процессов (подобно тому, как это делали авторы [8], но без переключений), устойчивость данных схем при моделировании решений SDEwMSs с разными начальными условиями, шагом дискретизации, рассматривая различные варианты задания матрицы переходов и меняя вре-
менной интервал. Будем использовать теорию и методы из [8], а также некоторые теоретические сведения из [1, 2]. Представленные в данной статье примеры могут быть полезны для понимания некоторых особенностей процесса моделирования систем с переключаемой диффузией.
2. Предварительные сведения
Пусть (О, F, Р) - вероятностное пространство;
Fh 10 < t < t0 + Т, - неубывающее семейство о-подалгебр F; ®г(-), г = 1, ..., d, - независимые винеровские процессы. Пусть М = {1, ..., т} - конечное множество.
Рассмотрим стохастическое дифференциальное уравнение с марковскими переключениями в форме
(1) dX(1) = а(Р(1),X(1))«1 + ]Т аг (3(1),X(1))«аг (t),
г=1
где ДО - однородный марковский процесс со счетным множеством состояний М;
(2) Рр(1 + Д) = /|/3(1) = и, х(5), 3(5), s < 1 )= qul(1 )Д + о(Д), и ФI,
где х(1) еЖп, а(-,-): Ж хМ ^ Ж и ст(у): Жп хМ ^ Жпхп. 2(-): Жп ^ Жтхт - ограниченная и непрерывная функция; Q(х) = (ди1 (х)) еЖтхт для каждого х; qul(x) > 0 при и Ф1,
т
X qul (х) = 0 для каждого и е М.
I=1
Предполагается, что функции Ав(1), х(1)) и стг(в(0, х(1)) определены при 1 е [10, 10 + Т], хеЖп , и удовлетворяют следующим условиям:
условию Липшица при всех 1 е [10, 10 + Т], х еЖп , у еЖп , и е М:
(3) |а(и, х) - а(и,у) + |ст(и,х) -и(и,у) < К|х - у|, а также условию
(4) |а(и, х)| + |ст(и,х)| < К(1 + |х|),
которое накладывает ограничения на скорость изменения компонент функций по х.
Будем использовать далее следующие обозначения:
|х| означает евклидову норму вектора х; ху - скалярное произведение векторов х и у; К - положительная константа.
Хих(0 или просто Х(0 - решение уравнения (1), удовлетворяющее начальным данным Хих(0 = х. Разобьем промежуток ^00, ^ + Т] точками деления к на N равных частей, так что tn+\ - tn = А, п = 0, 1, ..., N - 1. Приближение к Х(^) будем обозначать Уп, где Го = Х(^).
Представим схему Эйлера (Еиіег-Магцуата) в виде
(5) Гп+1 = Г + а(Рп, Г„ )А п + а(Рп, Г„ )А^, где р = ДО; х = х(0.
+ 2 °(рп , Гп )ст ' (рп , Гп )\.Аап )2 - Ап І где р = р(t), х = х(0.
Схему Тейлора (Платена) порядка 1,5 - в виде
(7) Гп+1 = Гп + а(Рп, Гп )А п + а(Рп, Гп )А^ +
+ 2 °(Рп , Гп )* '(Рп , Гп )|А®п )2 -А п }+
+ а ' (Рп, Гп МРп, Гп )А2п +
+ 2 (а(Рп, Гп )а '(Рп, Гп) + 2 а2(Рп, Гп )а "(Рп, Гп) ^ +
tп+1 ^2
гдеР = ДО; х = х(0; = | |dюSids2 ; Аак есть N(0, А) гаус-
Iп т
совски распределенных приращений винеровского процесса ш на подынтервалах tn < t <[п+\, = х0.
Используя леммы, доказанные Мильштейном, и обобщенную теорему Платена (см. [3]), можно установить среднеквадратичную сходимость схем (5), (6), (7) к решению уравнения (1) при выполнении условий (3), (4) соответственно порядка точности 1/2, 1, 3/2.
3. Теория постановки компьютерных экспериментов
Есть два основных типа задач, связанных с моделированием решений стохастических дифференциальных уравнений. Первый встречается в ситуациях, где требуется хорошее приближение, например в прямых моделированиях, фильтрации или тестировании статистических оценок. Во втором случае интерес сосредотачивается на приближении математических ожиданий функционалов процессов Ито, их вероятностных распределений и моментов [8].
Известно [2], что даже для численно устойчивой схемы при увеличении промежутка интегрирования ошибка численной схемы может расти и превышать допустимые пределы, хотя теоретически она остается ограниченной. Эта проблема становится особенно актуальной в тех задачах, где промежуток интегрирования ^0, Т] заранее не известен, например, при моделировании момента первого выхода. Однако в ряде случаев удается гарантировать выполнение требуемого ограничения сверху на ошибку численного метода. Как правило, достижение этого происходит за счет выбора подходящего шага интегрирования, параметров численной схемы и промежутка интегрирования. Существенную роль в устойчивости численных методов также играют особенности конкретных стохастических диффе-
ренциальных уравнении, к которым эти методы применяются
[1].
Пусть Yp, p = 0, 1, ..., - дискретная аппроксимация решения стохастического дифференциального уравнения Ито (1), которая стартует в момент t0 из точки Y0. Через Y’p, p = 0, 1, ..., обозначим эту же дискретную аппроксимацию, но стартующую в момент t0 из точки Y’0. Выберем некоторый конечный промежуток интегрирования [t0, T] и покроем его сеткой |т;. }nTq, такой
что
t0 = 1 <1 <... < TNT = T, л= max kj+1 -I\ .
T 0< j < Nt1 1
Для любого t e [t0, Т]определим натуральное число Nt следующим образом: Nt = max|/' -.ij <t}. [1]
Определение 1. Будем говорить, что дискретная аппроксимация Yp, p = 0, 1, ..., стохастически численно устойчива для данного стохастического дифференциального уравнения, если для любого конечного интервала [t0, T] существует такое постоянное число Д0 > 0, что Vs> 0 и VA e (0, Д0) выполняется условие [1]:
(8) IY -rUo SUP P К - YN\ - s }= 0.
|Y0-Y0|^0 t0<t<T
3.1. КРИТЕРИЙ АБСОЛЮТНОЙ ПОГРЕШНОСТИ
Критерий абсолютной погрешности есть математическое ожидание модуля разности между приближением и процессом Ито в момент T s = E(\XT - Y(T)|), что дает меру приближения на конце временного интервала [0, T].
Можно получить статистическую оценку абсолютной погрешности, используя компьютерные эксперименты.
Смоделируем N типовых траекторий процесса Ито и их приближений, соответствующих тем же самым типовым траекториям винеровского процесса, построенных с использованием той или иной численной схемы.
Обозначим как Хт,к и Yт,k значения к-х моделируемых траекторий в момент Т. Будем оценивать абсолютную погрешность как статистическую величину [8]:
1 N ,
(9) ё = -1 Хтк - ^ .
N к=1
3.2. ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ ДЛЯ АБСОЛЮТНОЙ ПОГРЕШНОСТИ
Для больших N, как известно из центральной предельной теоремы, погрешность е ведет себя асимптотически как гауссовская случайная величина и сходится по распределению к неслучайному математическому ожиданию е абсолютного значения погрешности при N ^ да.
Невозможно воспроизвести бесконечное число траекторий. Однако можно оценить среднее отклонение ст.? от е и затем использовать эти оценки при построении доверительного интервала для абсолютной погрешности е. Для этого нужно моделировать М партий по N значений в каждой и оценивать отклонение от е следующим образом. Обозначим как Yт,kj значение к-й воспроизведенной траектории в /-партии в момент времени Т и как ХТ£/ соответствующее значение процесса Ито. Усредненные погрешности
1 N | ,
(10) ё] = -1 ХтХ] - YтЛj\
N к=1 1
из М партий, ] = 1, 2, ...,М, являются тогда независимыми и приближенно гауссовскими для больших N.
Будем вычислять погрешности в каждой партии, потому что тогда мы можем использовать /-распределение Стьюдента для построения доверительных интервалов сумм независимых гауссовских случайных величин (или приближенно гауссовских) с неизвестной дисперсией.
В частности мы оцениваем средние величины для партий
1 М 1 М N
(11) е = — £ё, =-------££\Хтк1 -Ytk \
М ]=1 ] NM ]=1 ^ т,к,] а,]|
и затем используем формулу
чтобы оценивать дисперсию для найденных средних величин партий. Эксперименты показали, что средние значения могут интерпретироваться как являющиеся гауссовскими для партий размером N > 15; но обычно используют N = 100 [8].
100 (1 - а) % доверительный интервал имеет вид (г - Дг,£ + Де) с
где значение t1-aM-1 определяется из t-распределения Стьюдента с M— 1 степенями свободы. Для M = 20 и а = 0,1 мы имеем t1-aM-1 ~ 1,73 из таблицы (см. [8]). В этом случае абсолютная погрешность £ попадет в соответствующий доверительный интервал с вероятностью 1 - а = 0,9.
3.3. ЗАВИСИМОСТЬ АБСОЛЮТНОЙ ПОГРЕШНОСТИ ОТ РАЗМЕРА ШАГА
Размер шага, конечно, имеет влияние на величину абсолютной погрешности. Это видно более явно, если построить зависимость log2 £ от log2 Д. Известно, что функция ДД) = A Ду становится линейной в логарифмических координатах, и выполняется loga ДД) = loga A + у loga Д для логарифмов с основанием a Ф 1. В сравнительных исследованиях удобно выбирать шаг в форме Д = а- для k = 1, 2, ... и a > 1 [8].
4. Постановка задачи
Для практических моделирований определяем (см. [1, 2,
3, 8]) Дю = £л/Д , где £ — гауссовская переменная.
Перепишем схему (5):
(14) Yn+1 = Y„ + a(Pn, Yn)Д„ + a(Pn, Yn)4Цп .
В таком виде схема рассмотрена в [25, 27] и др., где представлено доказательство ее сходимости к решению уравнения dx(t) = а(Р(ґ), х(0) dt + о(Д0, х(0) dй>(t),
х(0) = х0; в(0) = во, где в - марковский процесс переключений. Перепишем (6) в виде
(15) Уп+і = ¥„ + а(Рп, Уп )А „ + а(Рп, Уп )1т +
+ а(Рп , Уп )ст ' (Р-> х) 1 (1,1) ,
где 1(1), 1(11) - соответствующие интегралы Ито.
Перепишем (7) в виде
(16) 7п+1 = Уп + а(Рп, Уп )А п + а(Рп, Уп )1т +
+ ^(Рп, Уп )а '(Рп, Уп) I а,1) + а '(Рп, Гп )^(Рп, Уп )1а,о) +
+ Га(Рп , Уп )а'(Рп , Уп ) + 2 ^ 2(Рп , Уп )а " (Рп , Уп ) V(0,0) +
+ (а(Рп, Уп )а" (Рп, Уп) + 2 °ЧРп, Уп )*" (Рп, Уп)) I (0,1) + + СТ(Рп , Уп )(<т (Рп , Уп )с " (Рп , Уп ) + (<Т ' (Рп , Уп ))2 )/(1,1,1),
где 1(1) = А®п; /(1,1) = 2 (А®п2 - А п); 1(о,о)= 1 А2п;
1 ____ ^ 1 Л
1 (1,0) = 2^А п ■ ^ + “^ ^2 ; где £ь £ - гаУссовские слУчай-
43
/(0,1) = А®пА п -1 (1,0); 1 (1,1,1) = 2 ^3 А®п2 -А п ^А®п ;(см. [3]).
ные величины;
1 [1
2 ^ 3
Далее представим результаты компьютерных экспериментов для нахождения решений SDEwMSs, используя численные схемы (14), (15), (16), и сравним качество аппроксимаций, используя изложенные выше практические методы.
5. Пример
Рассмотрим процесс Ито X = {X, / > 0}, удовлетворяющий линейному стохастическому уравнению
(17) dXt = f (Д, X,)Х^ + ^ (Д, X,)Х^ ,
на временном интервале [0, Т], X0 = 1. Пусть М = {1, 2, ..., га} -число состояний марковской цепи; Q - матрица интенсивностей; Р - матрица переходных вероятностей, Р = I + QД, где I - единичная матрица (см. [3, 27]); ДД, X,) и g(Д, X,) - непрерывные, ограниченные, дважды дифференцируемые функции, удовлетворяющие условиям (3), (4).
В численных схемах (14), (15), (16) а(Д, х) = ДД, х,) х, а(Д, х,) = g(Д, х,) х, для рассматриваемого уравнения
(( 1 Л >
V
(18) Xt = Xо ехр I у(Д, х) - - g (Д, х,) к + g(Д, х,К
является решением уравнения (17) для , е [0, Т ] и данного вине-
ровского процесса о = {а,, , > 0}.
Чтобы смоделировать линейно интерполированную траекторию аппроксимаций (14), (15), (16), зададим начальное значение Го = X0 и будем рекурсивно генерировать значения Уп с равным значением шага Д. Будем моделировать 20 партий по N = 100 траекторий в каждой.
Расчеты погрешности и ширины доверительного интервала будем проводить для разного размера шага дискретизации. Кроме того, учитывая переключаемый характер диффузионного процесса, экспериментально проследим, насколько устойчива та или иная схема к колебаниям значений коэффициентов в уравнении (17), изменению матрицы переходов на различных временных интервалах при разной интенсивности скачков.
Случай 1.
Пусть У(вь х) принимает два значения {аь а2}, соответствующие первому и второму состоянию марковской цепи. Пусть g(fiь х,) принимает два значения {Аь Л2}.
1. в =
; Р = I + вД;
а! = 20 + С08 х; а2 = 8т х + 5со8 2х; Х\ = 0,2х; Х2 = 0,08х.
при А = 0,001 (Т = 0,01); марковская цепь
БсїІаЬ СгарЫс (
Еііе ІРРІ5 ЫН:
£>1э1оМ ь|
. Ббігопд ТауІ« эсИете
Рис. 3. Аппроксимация Тейлора 1,5 и точное решение уравнения
2. Изменим значения а1, а2, Я1, Я2: а1 = соб2(х); а2 = Бт(х) + 5соб(2х); Л1 = 0,2х; Я2 = 0,5х.
0.0-1---1---1---1—---1---1---1---1—---1---1—-—
О 10 20 30 40 50 60 70 80 90 100
--- XI XI XI О б
— ;1 ^1 1 л1;1 і1А1;1;1 ь
А ,
Рис. 4. Аппроксимация Эйлера (красная кривая), Мильштейна (синяя кривая), Тейлора 1,5 (зеленая кравая) и точное решение уравнения А = 0,025 (Т = 1)
Рис. 5. Среднее значение абсолютной ошибки и доверительный интервал при Д = 0,2; 0,1; 0,05; 0,025 (Т = 1); зависимость погрешности от размера шага дискретизации
Scilab Graphic (0)
File Tools Edit
£>|э|ок»Ы
ytl
Рис. 6. Аппроксимации Эйлера, Мильштейна, Тейлора 1,5 и точное решение уравнения Д = 0,00375; (Т = 0,3)
Рис. 7. Среднее значение абсолютной ошибки и доверительный интервал при Д = 0,03; 0,015; 0,0075;0,00375 (Т = 0,3)
ЗсНаЬ СгарЫс (0)
11°) |^3д1
Н1е Хоок Е(№
Д1|э|о|сер| Ь|
. 5 з&опд Таукх эсЬегие
Рис. 8. Аппроксимации Эйлера, Мильштейна, Тейлора 1,5 и точное решение уравнения Д = 0,0005; (Т = 0,02)
Рис. 9. Среднее значение абсолютной ошибки и доверительный интервал при Д = 0,004; 0,002; 0,001; 0,0005 (Т = 0,02)
Случай 2.
1. б
ґ- 5 sin х 0 0
7cos х
- cos2 х
0
5sin х 0
- cos 2х cos 2 х
0
cos2 х
7 cos х
а1 = 3sin х; а2 = 2 + cos х; а3 = -6sin 2х; а4 = 0,2 х; Х1 = 0,2х; Х2 = 0,07х; Х3 = 2х2; Х4 = 0,8х.
0
Изменим значения:
а1 = sin x; а2 = 1 + cos x; а3 = 2sin 2x; а4 = 0,2x2; X\ = 0,2x; Я2 = 0,7x; X3 = 0,1x2; A4 = 0,5x;
точное решение уравнения при А = 0,025; (T = 1)
интервал при Д = 0,2; 0,1; 0,05; 0,025 (Т = 1); зависимость погрешности от размера шага дискретизации
Рис. 14. Среднее значение абсолютной ошибки и доверительный интервал при Д = 0,03; 0,015; 0,0075; 0,00375 (Т = 0,3); зависимость погрешности от размера шага дискретизации
интервал при А = 10-1; 10-2; 10-3; 10-4 (Т = 1; 0,1; 0,01; 0,001)
2. б =
- 5бш х
0
2 2 СОБ х - СОБ х
5бш х 0
0 - соб2 х соб2 х
0
0 0 0
7соб х - 7соб х^
а1 = 2 + біп х; а2 = 1 + соб х; а3 = 3 соб2 2х; а4 = біп 4х; Х1 = 0,5х; Х2 = 3х; Х3 = 0,2х; Х4 = 0,007х.
0
File Tools Edit
Э|Э|о1«а| k\
Рис. 17. А = 0,025; (T = 1)
Scïlab Graphie (0) File look Edit
ddJiüLi"
Рис. 18. A = 0,1; 0,05; 0,025 (T = 1)
Рис. 19. Аппроксимации Эйлера, Мильштейна, Тейлора и точное решение уравнения при Д = 0,01; (Т = 0,02)
Рис. 20. Среднее значение абсолютной ошибки и половина ширины доверительного интервала при А = 0,02; 0,01 (Т = 0,2); зависимость погрешности от шага дискретизации
Далее используем Л1 = 0,05; Я2 = 8; Я3 = 0,2; Я4 = 0,007.
epsl =
96.423884 2.1340132 0.2344813 0.0532120 differl =
eps2 =
5442.6701 0.8128139 0.1081653 0.0114317 differ2 =
eps =
Nan 0.7232948 0.0828786 0.0072688 differ =
27.83168 3335.0331 Nan
0.6632425 0.2606815 0.1054961
0.0542171 0.0193577 0.0152877
0.0152258 0.0021111 0.0013275
Рис. 21. Среднее значение абсолютной ошибки и половина длины конфиденциального интервала при А = 0,1; 0,02; 0,004; 0,0008 (T = 1; 0,2; 0,04; 0,008); средняя погрешность в партии при А = 0,0008 (T = 0,008); марковская цепь при А = 0,0008
(T = 0,008)
а1 = sin x; а2 = -4 cos x; а3 = 3 cos2 x; а4 = cos 2x;
Ai = 0,2; Я2 = (3 + 0,1x); A3 = 0,7; Я4 = 0,8.
Scilab Graphic (0)
Рис. 22. Аппроксимация Мильштейна при А = 0,02 (T = 0,1)
Єр5І =
0,2184561 0,2083285 differl -
0.0209476
0,0190861
ерь2 =
0.1792578 0.1612028 ¿і ГҐЄГ2 =
0.0133727
0.0133734
ерз =
0.1719160 0.1610206 differ =
0.0136489
0.0179922
Рис. 23. Среднее значение абсолютной ошибки и половина длины доверительного интервала при А = 0,01; 0,02; величина средней погрешности в партии А = 0,02 (Т = 0,1)
Случай 3.
( 3
б =
cos X 0 0 0 0 - 3 С08 х 0 0 "
0 С08 х 0 - С08 х 0 0 0 0
0 0 0- - 5 С08 2 х 5 С08 2 х 0 0 0
0 0 0 - 6 8ІП 2х 0 0 0 6 8Іп 2 х
0 0 0 - 2 8ІП х 0 28ІП х 0
0 - 5 С08 х 0 0 5 С08 х 0 0
2 8Іп 4 х 0 0 -28Іп4 х 0 0 0
0 СЧ х 3 - 0 0 0 2 х 3 0 ,
0 0 0 Ч0
а1 = 3х; а2 = 2sin х; а3 = -6х2; а4 = cos 4х; а6 = ^Іп х; а7 = 2 + х; а8 = 5sin 2х;
Х1 = 0,1х; Х2 = 0,03х; Х3 = 7х; Х4 = 0,7х;
Х5 = 0,9х; Х6 = 0,2х; Х7 = 4х; Х8 = 0,2х.
а5 = cos х;
точное решение уравнения А = 0,00375 (Т = 0,3); марковская цепь
А1 = 0,1х; А2 = 0,3х; А3 = 0,7х; А4 = 0,7х; А1 = 0,9х; А2 = 0,2х; Я3 = 0,4х; Я4 = 0,2х.
£ІІе ТооЬ ЕА
I | О I 6Ев| ^иіег зс(іегг»е
Рис. 26. Среднее значение абсолютной ошибки и доверительный интервал при А = 0,05; 0,025; 0,0125; 0,00625 (Т = 0,5)
б =
3соб х 0 0 - 3 соб х 0 0 0 0
0 соб х 0 0 - соб х 0 0 0
0 0 0 5соб2 х 2 - 5соб х 0 0 0
6біп2 х 0 0 0 0 0 0 - 6 БІП 2 х
0 0
0 0
0 2 біп 4 х 00
0
- 5 соб х 0 0
0
0
0
- 3біп х2
- 2біп х
0
- 2 біп 4 х
0
0
5соб х 0 0
2біп х 0 0 0
0
0
0
3біп х 2
а1 = 3х; а2 = 2 біп 4х; а3 = 6 соб х; а4 = біп 4х; а5 = соб х; аб = -біп х; а7 = 2 + х; ав = 5 біп 2х;
Х1 = 1,5х; Х2 = 3х; Х3 = 0,7х; Х4 = 4х;
Х5 = 0,29х; Х6 = 0,6х; Х7 = 0,24х; Хв = 0,27х.
точное решение уравнения А = 0,0025(Т = 0,1)
Рис. 28. А = 0,01; 0,005; 0,0025 (Т = 0,1)
точное решение уравнения А = 0,00625 (Т = 0,5)
зации (А = 0,05; 0,025; 0,0125 (Т = 0,5))
Рис. 32. Аппроксимации Мильштейна (синяя кривая), Тейлора (зеленая) А = 0,025 (T = 1)
Рис. 33. Среднее значение абсолютной ошибки и доверительный интервал схем Эйлера, Мильштейна (на графике), Тейлора (в окне программы) при А = 0,2;0,1;0,05; 0,025 (Г=1)
6. Заключение
В данной статье были использованы известные численные схемы Эйлера, Мильштейна, Тейлора (Платена) в виде (14),
(15), (16) для нахождения приближенного решения SDEwMSs вида (17) с механизмом переключений, зависящим от фазового состояния.
Используя методы Peter E. Kloeden, Eckhard Platen, Henri Schurz из [8], в данной работе было проведено сравнение качества аппроксимаций, подсчитаны значения абсолютной
погрешности, ширины 90% доверительного интервала с различными вариантами выбора начальных условий, шага дискретизации, ширины временного интервала, размера и вида матрицы переходов.
В результате можно сделать следующие выводы.
Качество аппроксимаций зависит от вида функций ДД, Х) и g(Д, Х), являющихся коэффициентами в уравнении (17) и для которых требуется выполнение условий (3), (4). Эксперименты показали, что данные условия можно несколько ослабить, но в таком случае качественные приближения возможны только на малых временных интервалах с малым размером шага дискретизации.
В целом точность методов Тейлора и Мильштейна определенно выше на небольших временных интервалах с малым размером шага дискретизации. При увеличении временного интервала разница в ошибке уменьшается и качество аппроксимаций мало отличается.
Устойчивость все три метода показывают на малых временных интервалах. На достаточно больших временных интервалах часто ошибка численной схемы превышает допустимые значения. В любом случае, чрезмерное уменьшение шага дискретизации приводит к увеличению погрешности. Накопление ошибки происходит из-за того, что на каждом шаге винеровский процесс получает некоторое приращение, и, соответственно, при малом размере шага и относительно большой длине временного интервала приращение Аш становится настолько большим, что не компенсируется малым размером шага А, как предполагалось в теоретических исследованиях. Схема Эйлера несколько более устойчива в таких случаях. Схема Мильштейна, отличаясь от схемы Эйлера только одним слагаемым, часто показывает более точные результаты по сравнению со схемой Эйлера, и в тоже время она немного устойчивее «чувствительной» к увеличению интенсивности и диапазона значений скачков схемы Тейлора (Платена).
Интересно было проследить, как меняется качество аппроксимаций при увеличении размера матрицы переходов. Это
также одно из возможных направлений для теоретических исследований.
Отметим также, что при проведении экспериментов по моделированию сложных процессов многое зависит от конкретных уравнений, начальных условий, корректности решения (в нашем случае само «точное решение» тоже находится несколько приближенно), выбора интервала, шага дискретизации. Поэтому однозначных выводов о качестве методов по нескольким примерам сделать нельзя. Показанные в данной статье иллюстрации численного моделирования представляют лишь небольшую часть реально проведенных компьютерных экспериментов и демонстрируют часто повторяющиеся случаи, которые позволяют все же сделать некоторые выводы о моделировании подобных процессов и могут помочь при решении практических задач.
Литература
1. КУЗНЕЦОВ Д.Ф. Стохастические дифференциальные уравнения: теория и практика численного решения. - Спб.: Издательство Политехнического университета, 2007. -800 с.
2. МИЛЬШТЕИН Г.Н. Численное интегрирование стохасти-
ческих дифференциальных уравнений. - Свердловск: Изд-во Уральского университета, 1988. - 224 с.
3. ЧЕРНЫХ Н.В., ПАКШИН П.В. Алгоритмы численного решения стохастических дифференциальных систем с переключаемой диффузией // Управление большими системами. Сборник трудов. - 2012 - №36 - С. 106-143.
4. BAO J., HOU Z. An analytic approximation of solutions of stochastic differential delay equations with Markovian switching // Mathematical and Computer Modelling. - 2009. - Vol. 50, №9-10. - P. 1379-1384.
5. DARGATZ C. A Diffusion Approximation for an Epidemic Model // Ludwig-Maximilian University Munich. - 2007. -
P. 1- 20. - URL: http://www.epub.ub.uni-muenchen.de/1882/ (дата обращения: 19.11.2012).
6. ELLIOTT R.J., KRISHNAMURTHY V., SASS J. Moment based regression algorithms for drift and volatility estimation in continuous-time Markov switching models // Econometrics Journal - 2008. - Vol. 11. - P. 244-270.
7. FAN J., ZHANG C. A Reexamination of Diffusion Estimators
With Applications to Financial Model Validation // Journal of the American Statistical Association. - 2003. - Vol. 98, №461. - P. 118-134.
8. KLOEDEN P.E., PLATEN E., SCHURZ H. Numerical Solution of SDE Through Computer Experiments. - Berlin: Springer -Verlag, 1994. - 294 p.
9. LI H., XIAO L., YE J. Strong predictor-corrector Euler-
Maruyama methods for stochastic differential equations with Markovian switching // Journal of Computational and Applied Mathematics. Available online 5 July 2012. - URL:
http://www.pas.sciplore.org/archivepages/2011-03/ (дата обращения: 19.11.2012).
10. LI R., PANG W.-K., LEUNG P.-K. Convergence of numerical solutions to stochastic age-structured population equations with diffusions and Markovian switching // Applied Mathematics and Computation. - 2010. - Vol. 216, №3. - P. 744-752.
11. LIU R.H., ZHANG Q., YIN G. Option pricing in a regime -switching model using the fast fourier transform // Econometrics Journal. - 2008. - Vol. 11. - P. 244-270.
12. MAO X., YUAN C. Stochastic Differential Equations with Markovian Switching. - London: Imperial College Press, 2006. -428 p.
13. MAO X., YUAN C., YIN G. Numerical method for stationary distribution of stochastic differential equations with Markovian switching // Journal of Computational and Applied Mathematics. -2005. - №174. - P. 1-27.
14. MAO X., YUAN C., YIN G. Approximations of Euler-Maruyama type for stochastic differential equations with Markovian switching, under non-Lipschitz conditions // Journal of
Computational and Applied Mathematics. - 2007 - Vol. 205, №2. - P. 936-948.
15. MARITON M. Jump Linear Systems in Automatic Control. -Marcel Dekker, New York, 1990. - 297 p.
16. MILOSEVIC M., JOVANOVIC M. A Taylor polynomial approach in approximations of solution to pantograph stochastic differential equations with Markovian switching // Mathematical and Computer Modelling. - 2011. - Vol. 53, №1-2. -P.280-293.
17. PREIS T., SCHNEIDER J.J., STANLEY H E. Switching processes in financial markets // Proceedings of the National Academy of Sciences (USA). - 2011. - Vol. 108, №19. -P. 7674-7678.
18. RATHINASAMY A. Split-step 9-methods for stochastic age-dependent population equations with Markovian switching // Nonlinear Analysis: Real World Applications. - 2012. -Vol. 13, №3 - P. 1334-1345.
19. RATHINASAMY A., YIN B., YASODHA B.. Numerical analysis for stochastic age-dependent population equations with Poisson jump and phase semi-Markovian switching // Communications in Nonlinear Science and Numerical Simulation. -2011. - Vol. 16, №1. - P. 350-362.
20. RILEY D, RILEY K. Simulation of stochastic hybrid system with switching and reflecting boundaries [Электронный ресурс] // IEEE Proceedings of the 2008 Winter Simulation Conference [Internet-Conference]. - 2008. - P. 804-812.
21. SIU T.K., LAU J.W., YANG H. Pricing Participating Products under a Generalized Jump-Diffusion Model // Journal of Applied Mathematics and Stochastic Analysis. - 2008. - P. 1-30.
22. SONG Q.S., YIN G., ZHANG Z. Numerical Solutions of Stochastic Control Problems for Regime-switching Systems // Submitted to Dynamics of Continuous, Discrete and Impulsive Systems. - URL: http://www.math.wayne.edu/~zzhang/ publist.html (дата обращения: 18.11.12).
23. WU S.J., ZHOU B. Existence and uniqueness of stochastic
differential equations with random impulses and Markovian
switching under non-lipschitz conditions // Acta Mathematica Sinica. - 2011. - Vol. 27, №3. - P. 519-536.
24. YIN B., MA Z. Convergence of the semi-implicit Euler method for neutral stochastic delay differential equations with phase semi-Markovian switching // Applied Mathematical Modelling. -2011. - Vol. 35, №5. - P. 2094-2109.
25. YIN G.G., ZHU C. Stochastic modeling and applied probability. Hybrid switching diffusions. Properties and applications. -Springer Science + Business Media, LLC, 2010. - 395 p.
26. YIN G.G., KRISHNAMURTHY V. LMS Algorithms for Track-
ing Slow Markov Chains With Applications to Hidden Markov Estimation and Adaptive Multiuser Detection // IEEE transactions or information theory. - 2005. - Vol. 51, №7. -P. 2475-2490.
27. YIN G., MAO X., YUAN C., CAO D. Approximation methods for hybrid diffusion systems with state-dependent switching processes: numerical algorithms and existence and uniqueness of solutions // SIAM Journal on Mathematical Analysis. -2010. - Vol. 41, №6. - P. 2335-2352.
28. YUAN C., MAO X. Convergence of the Euler-Maruyama
method for stochastic differential equations with Markovian switching // Mathematics and Computers in Simulation. -2004. - Vol. 64, №2. - P. 223-235.
MODELLING SOLUTIONS OF STOCHASTIC DIFFERENTIAL EQUATION WITH MARKOVIAN SWITCHINGS
Nadezda Chernykh, Arzamas Polytechnic Institute of R.E. Alekseev, Nizhny Novgorod State Technical University, Arzamas, post-graduate student ([email protected]).
Abstract: We consider mathematical models of hybrid systems in the form of stochastic differential equations with Markovian switchings and a state-dependent switching component. We propose an extension of Euler, Milshtein, and Taylor schemes for numerical approximation of solutions for such stochastic differential equations. Quality of approximations is experimentally compared in Scilab software for numerical computation.
Keywords: stochastic systems, Markovian switchings, Euler scheme, Milstein scheme, Taylor scheme, convergence, stability, error, confidential interval.
Статья представлена к публикации членом редакционной коллегии Я. И. Квинто
