CC BY
203
ИСПОЛЬЗОВАНИЕ ПОСЛЕДОВАТЕЛЬНЫХ МЕТОДОВ МОНТЕ-КАРЛО ДЛЯ ОЦЕНИВАНИЯ РИСКОВ НА ФИНАНСОВЫХ РЫНКАХ
Наливкин Д. В.
(Липецкий государственный технический университет,
Липецк)
В статье предложена методика стохастического моделирования временных рядов финансовых рынков на основе модифицированного процесса Орнштейна-Уленбека. Предложен метод идентификации модели, опирающийся на последовательные методы Монте-Карло и принцип максимального правдоподобия. Рассмотрены практические результаты применения данной методики для вычисления значений Value-At-Risk с целью оценки рисков на рынке акций.
Ключевые слова: финансовый рынок, Value-at-Risk, процесс Орнштейна-Уленбека, последовательные методы Монте-Карло.
Введение
При исследовании поведения рынков в краткосрочном периоде актуальным является исследование колебаний (волатильности) финансовых рынков, так как в таких масштабах рынки демонстрируют существенную нестабильность. Исследователями разработан ряд моделей, которые позволяют моделировать поведение рынков как случайный процесс. Перечислим наиболее известные типы таких моделей:
1) Параметрические методы.
2) Модели условной скедастичности ARCH/GARCH.
3) Модели стохастической волатильности [1, 8].
Каждый из этих типов обладает определенными достоинствами и недостатками. Недостатками моделей первого и второго типа является плохая способность адекватно отражать резкие
колебания цен ввиду того, что на реальных рынках распределения цен не являются нормальными и имеют так называемые «тяжелые хвосты». Недостатком моделей третьего типа является излишняя трудность идентификации параметров модели, которая следует из того, что для таких моделей нельзя в явном виде выписать функцию правдоподобия.
Адекватность каждой модели определяется во многом тем, какие методы могут быть использованы для идентификации параметров модели и какой результат они дают. Как правило, для идентификации сложных моделей применяются методы Монте-Карло. Однако большинство этих методов, являются малоэффективными в вычислительном плане и требуют больших вычислительных затрат.
С другой стороны, в последние годы исследователями интенсивно разрабатывается класс методов Монте-Карло, основанных на принципах байесовского анализа, которые в определенном круге задач более эффективны по вычислительным затратам и дают лучшие результаты, чем стандартные подходы. Эти методы называются последовательными методами Монте-Карло или методами частичного фильтра (particle filter) [5]. Соответственно, имеет смысл их применение для моделирования и идентификации моделей стохастической волатильности, что было сделано некоторыми исследователями. Тем не менее, ими не рассматривался вопрос построения таких моделей, которые позволили бы максимально использовать потенциал последовательных методов Монте-Карло за счет соответствующим образом спроектированной структуры.
В данной статье предложена стохастическая модель движения цен на рынке на основе процесса Орнштейна-Уленбека, в которой параметры процесса представляют собой случайные величины с неизвестным заранее нестационарным распределением. Такая модель весьма эффективно реализуется через модификацию метода APF (auxiliary particle filter, последовательного частичного фильтра) [7] из класса последовательных методов Монте-Карло одновременно с последовательной идентификацией распределений параметров модели. Идентификация параметров опирается на аналитические формулы для оценки па-
раметров, полученные из принципа максимального правдоподобия для процесса Орнштейна-Уленбека. Практика показала, что подобный подход, с одной стороны, достаточно экономичен в вычислительном плане, с другой стороны, позволяет строить модель распределения рыночных цен, которая адекватно отражает «тяжелые хвосты» распределений цен и адаптируется под резкие колебания цен. Анализ значений Value-at-Risk, вычисленных по результатам исторического моделирования для различных акций российского рынка, подтвердил практическую значимость данной модели.
1. Постановка задачи
В качестве базовой модели рассматривается модель, являющаяся модификацией средневозвратного стохастического процесса Орнштейна-Уленбека. Стандартный вид уравнения этого процесса c ненулевым средним выглядит следующим образом [1, 3, 6]:
(1) dx = Х(р - x)dt + odW,
где dW — броуновское приращение, f, X, a — параметры модели (среднее, коэффициент дрейфа и волатильность соответственно). Исследуемая в данной работе модель отличается от стандартной тем, что параметры модели f, X, а, рассматриваются как случайные величины, независимые от x(t) и распределенные по заранее неизвестному нестационарному закону. Принимая, что логарифмы значений цен log(y(t)) являются реализацией случайного процесса x(t), измеренного в последовательные моменты времени с некоторым постоянным интервалом At, получаем общий вид модели движения цен, рассматриваемой в данной работе:
dx = A(t)(/u(t) - x(t))dt + o(t)dW, log(j(t)) к x(t)
( ) A(t) к Л^), /u(t) к M(t),o(t) к E(t)
Выбор процесса Орнштейна-Уленбека в качестве основы для предлагаемой модели обусловлен тем, что данный процесс, является наиболее простым из класса средневозвратных процессом и поэтому часто применяется в исследованиях финансовых рынков, которые демонстрируют аналогичное поведение в крат-
косрочной перспективе. При этом он имеет аналитическое решение задачи определения точечных оценок параметров по реализации процесса. Это решение опирается на нахождение минимума функции максимального правдоподобия для аналитического решения стохастического уравнения (1) [3]. Если даны значения реализаций процесса х1, х2, хп в моменты времени / = 0, Д/, 2Д/, пД/, то значения этих оценок ц , А , <г могут быть вычислены по следующим формулам:
]И =
- Я у) - (Я 2 - ^ ^у)
х £.
/ =-----------1п
1 , Яху - №х - №у + ПИ2
2
(3) А «хх - 2 &х + п£
а2 =----2/1 . [«уу - 2бсБ1у + а2-
п -2ЛА^ уу ху хх
п(1 - е )
-2^(1 -а)(§у -а&х) + п^2(1 -с?)2],с? = е~АА‘ где величины « представляют собой набор статистик:
= &-„ «у =±хг, =]Гхг^1,
(4) г=1 г=1 г=1
^ху =]^хг-1 х , «уу =^их1
1=1 1=1
Вычисления по этим формулам являются наиболее простым способом получения точечных оценок параметров процесса. Тем не менее, если параметры являются нестационарными, то прямое использование вышеприведенных формул для расчетов невозможно, так как они дают единичные значения. Наиболее простым решением данной проблемы является последовательное использование для идентификации параметров в моменты времени к = 0, 1, ..., п по реализациям процесса за предыдущие N моментов времени: к - 1, к - 2, ..., к - N где N - константа. Однако этот метод чреват существенными проблемами с выбором числа N - при излишне больших значения идентифицированные параметры будут сглаженными и невозможно будет 74
корректно идентифицировать резкие колебания параметров, при излишне малом N точность вычислений по формулам (3)-(4) резко падает и результат становится неприемлемым. Наконец, вычисления по формулам (3)-(4) дают точечные оценки параметров, что недостаточно для идентификации распределений параметров в предложенной автором модели. Тем не менее, эти соотношения могут быть использованы для последовательной идентификации распределений параметров в рамках методов Монте-Карло.
2. Моделирование
Классической методикой, используемой повсеместно для решения задач моделирования и идентификации стохастических процессов, является метод марковских цепей Монте-Карло. Его суть заключается в построении марковского процесса, распределение которого сходится к распределению стохастического процесса. Недостатком этой методики является низкая вычислительная эффективность и невозможность определить заранее число шагов моделирования, необходимое для сходимости метода. Для преодоления этих недостатков в последние десятилетия был разработан ряд методов, опирающихся на применение методик байесовского анализа для повышения эффективности и получения более точных результатов в соответствующих задачах. Если подобные методы правильно спроектированы, то они могут давать гораздо лучший результат, чем стандартные методы Монте-Карло. Кроме того, они также позволяют применять байесовские методики для идентификации параметров моделей или их распределений, что также является плюсом по сравнению с классическими методами, которые, как правило, ограничены использованием поисковых или итеративных методов для получения точечных оценок параметров. Одним из наиболее известных методов данного класса является метод «частичного фильтра» (particle filter) [5], известный также как последовательный метод Монте-Карло. Этот метод можно описать следующим образом:
Пусть имеется последовательность наблюдаемых значений yk в моменты времени k = 0, 1, 2, ..., T, которые являются реали-
зациями неизвестного распределения p(xk). Задачей является построение последовательности распределений p(xk) в соответствующие моменты времени. Приближение xk, которое может быть построено с использованием байесовских методик, будет иметь вид апостериорного распределения p(xk | y0, y1, ., yk).
Методы частичного фильтра опираются на предположение, что xk и yk могут быть смоделированы в следующей форме:
1) Значения x0, x1, . , xk являются марковским процессом первого порядка:
(5) xk \ xk-1 K Pxk|xk-1(x I xk-1)
с начальным распределением p(x0);
2) Наблюдаемые значения y0, y1, ., yk условно независимы при условии, что x0, x1, ..., xk известны. Другими словами, каждое yk зависит только от xk:
(6) yk 1 xk х py\x (y1 xk).
Методы частичного фильтра, так же как и все методы Монте-Карло, опираются на создание множества реализаций случайного процесса. В данном случае множество реализаций аппроксимирует распределение p(xk | y0.k), которое представляется в дискретном виде как множество точек (частиц) с весами:
(7) {wkL), xkL)): L = 1,...,р}.
Веса значимости wkL) являются приближениями относительных апостериорных вероятностей частиц, соответственно, их сумма равна единице.
За основу в данной работе был взят алгоритм APF (auxiliary particle filter, последовательного частичного фильтра) из [7]. Этот алгоритм может быть описан следующим образом.
Пусть задано некоторое начальное распределение xk в нулевой момент времени p(x0).
Для k = 1, 2, ., T последовательно выполняются следующие шаги:
1) Для i = 1, ..., Р вычисляется в(()— значение, ассоциированное с распределением p(xk+1 | x^), например, математическое ожидание, из которого берется частица с индексом i;
2) Для і = 1, Р вычисляются веса первой ста-
= **Й
дииИк(1) = р(ук | в(()) и нормализуются:
1)
і=1 •
?
3) Производится повторная выборка с целью получения нового массива \х^, <; : Р раз моделируется дискретная
случайная величина со значениями £ = 1, ..., Р и вероятностями значений Ик(г), к = 1, ..., Р. Частицы, которые имеют индексы полученных реализаций £, формируют новый
массив частиц {4і ^ Г=,;
4) Для і = 1, ..., Р производится распространение частиц через выборку из плотности вероятностей переходов: хк) к Р(хкг) I хк-11,^г), и вычисляются истинные веса зна-
(і )>
чимости: <> .
‘ мл №">>
Вышеизложенный алгоритм справляется с задачей моделирования случайного процесса, но в нем не заложена возможность идентификации параметров. Однако данная модель позволяет разработать расширение данной вычислительной схемы с последовательной идентификацией параметров.
Автором предложена модификация алгоритма АРБ с после -довательной идентификацией параметров для процесса Орн-штейна-Уленбека, которая может быть построена следующим образом.
Рассмотрим в качестве базовых элементов модели частичного фильтра не единичные значения {хк}, а расширенный вектор {хк, ак, Хк, рк}. На каждом этапе работы частичного фильтра будем обновлять значения ак, Хк, рк по формулам (3), используя для этого набор статистик {^х, Sx, ^хх, ^ху, 5^}. Эти статистики строятся по N предыдущим значениям хк_1, хк-2, ..., хк_м, где N -заранее определенная константа. Соответственно, одновременно с аппроксимацией распределения х можно аппроксимировать
апостериорные распределения параметров с использованием аналогичных соотношений:
|Ч">, хк">): Ь = 1,..., р} М(Г > 4Ч">, /ь): Ь = 1,..., р}
(8) 1(0 *{ЧЬ>,^Ь>):Ь = 1,...,р}
Л(? >*{>,4 )ь = 1,..., р}
Далее, распространение частиц через выборку из плотности переходов в шаге 4 алгоритма опирается на моделирование процесса Орнштейна-Уленбека в интервале времени от к до к + 1. Благодаря тому, что процесс Орнштейна-Уленбека протекает в непрерывном времени, в заданном интервале процесс можно моделировать в виде последовательности реализаций с неким временным шагом 5 = ИЯ, где Я - выбранное целое неотрицательное число. Соответственно, процесс распространения частиц с момента времени к до момента времени к + 1 можно аппроксимировать как результат последовательного применения Я
шагов следующего вида к исходному набору частиц {хк’> }=1:
х£+1/Я = Я4>(/к> - хк°> + >N(0,1>
(9> хк+;/Я = Я4>(/4° -х£(;-1>/Я> + ак~>N(0,1>
х& = Я4>(/« -х«(Я_1>/я> + °к>N(0,1>
где N(0, 1> - реализация стандартного нормального распределения. Выбор Я оказывает существенное влияние на результат работы алгоритма.
Таким образом, автором разработан следующий алгоритм, который представляет собой модифицированную версию алгоритма АРБ для процесса Орнштейна-Уленбека вида (2> с байесовским приближением распределений параметров модели.
Пусть заданы начальное распределение процесса р(х0) и начальные значения распределений параметров р(оь>, р(и0>, р(Х0>.
Для к = 1, 2, ..., Т последовательно выполняются следующие шаги:
1> Для 1 = 1, ., Р вычислить в£> = 4к2>(/к> - х!( >>.
2> Для 1 = 1, ..., Р вычислить веса первой стадии: ^> = ^->1 Р(Ук I вк>>, и нормализовать их:
Ж(’>
И?" = -р—
1 ’
1=1 .
3> Производится повторная выборка с целью получения нового массива {хк(г>,%г }=1 : Р раз моделируется дискретная случайная величина со значениями £ = 1, ..., Р и вероятностями значений > , к = 1, ..., Р. Частицы, которые имеют индексы полученных реализаций £ , формируют новый массив частиц {хк(1 >,^}=1 . Соответственно, также формируются наборы значений параметров модели
{,4 >, /'>}
4> Для 1 = 1, ..., Р производится обновление значений частиц {хк+1 }=1 по формулам (9>, получение истинных весов зна-
Р( Ук 1 ~к(1>> к Р( У к |вк^' >>’
чимости по формуле: Жк(1> = -
5> Обновление значений статистик через формулы (4> по новым значениям {х(1_+1 }=1. Вычисление новых значений параметров {тк^, 4+1, /к1+1}=1 на основе новых значений статистик по формулам (3>.
Используя этот алгоритм, можно построить на основе ряда исторических данных У0, У1, . , УТ аппроксимированное распределение логарифмов цен р(хк\ у0:к> для каждого момента времени к = 1, 2, ..., Т. Эти распределения можно использовать для вычисления мер риска Уа1ие-а1>Ш8к в некотором временном горизонте (к + 1, ..., к + 0. Они могут быть вычислены как квантили распределений р(хк+11 У0.к>, •••,р(хк+\ У0.к>. Для этого необходимо произвести моделирование случайного процесса в промежутке времени от к до к + t по формулам (9>. Но параметры процесса
изменяются во времени, соответственно, если принять, что они будут неизменны на отрезке моделирования (k + 1, ..., k + t), то точность вычисленного значения Value-at-Risk будет весьма низка. Наиболее простым решением данной проблемы является
линейная интерполяция значений параметров й4,^'), ft<,) L:
(1т +вЛ=,М1 15Г.1={,)+/vL •
( ) ki
и подстановка получающихся значений в формулы (9) при моделировании процесса на интервале от k до k + t. Значения коэффициентов Дт, вц, вь вычисляются по известным формулам линейной интерполяции на основе предыдущих значений параметров в моменты времени k, k- 1, ..., k- C, где C - некоторая константа. Вычислительные эксперименты подтвердили обоснованность данного подхода.
3. Результаты
В качестве вычислительных экспериментов для оценки качества полученной модели было произведено историческое тестирование значений Value-At-Risk для данных по дневным котировкам акций российских компаний (таблица 1, таблица 2). Величиной, которая характеризует адекватность значений Value-At-Risk, является фактически наблюдаемый доверительный интервал:
(11) анабл = Р(У > Var (a, t) = NKopp / N ,
где N - общее количество интервалов, для которых был вычислен Value-at-Risk; NKopp - количество интервалов расчета Value-at-Risk, в которых максимальное падение цены было больше рассчитанного значения; а и t -заданный уровень доверия и временной горизонт соответственно. Чем меньше анабл отклоняется от заданного уровня доверия, тем адекватнее оценка рисков. В экспериментах были использованы следующие параметры: количество частиц P = 2000, количество промежуточных шагов R = 4, размер выборки для идентификации параметров
N = 7, размер выборки для линейного интерполирования параметров С = 7.
Таблица 1. Результаты исторического тестирования для котировок акций РАО ЕЭС за 2006-2007 гг._____________________
Уровень доверия
Временной горизонт, дней 0,9 0,95 0,99
1 0,894 0,944 0,995
3 0,877 0,954 0,982
5 0,884 0,952 0,976
10 0,886 0,954 0,976
Таблица 2. Результаты исторического тестирования для котировок акций АО ГАЗПРОМ за 2006-2007 гг.__________________
Уровень доверия
Временной горизонт, дней 0,9 0,95 0,99
1 0,890 0,952 0,984
3 0,888 0,929 0,982
5 0,901 0,950 0,991
10 0,907 0,950 0,995
Следует сказать несколько слов о погрешностях вышеприведенных цифр. Объем выборки значений цен для вышеизложенного исторического тестирования составил во всех случаях 497 значений. В случае последовательных методов Монте-Карло моделирование производится последовательно во времени с количеством шагов, по времени равным размеру выборки (в данном случае 497), причем на каждом шаге Р раз производится выборка из приближаемого распределения, где Р - это количество "частиц" в алгоритме, в данном случае 2000. Соответственно, общее количество испытаний метода Монте-Карло можно
оценить в 497 • 2000 —106, что дает погрешность вышеприведенных цифр порядка 0,001. С другой стороны, погрешность можно оценить по приближенным формулам специально для метода "частичного фильтра" [4]. Эти формулы для данного случая дают величину порядка 2/2000=0,001, что сходится с первой оценкой. Эти значения меньше, чем погрешность, связанная с размером выборки для тестирования и составляющая ±1/497 — ±0,002, соответственно, именно это число (0,002) можно считать погрешностью вышеприведенных оценок анабл.
Как видно по данным, приведенным в таблице 1 и таблице 2, разработанный подход дает весьма точную оценку значений Value-at-Risk для различных акций, причем для уровня доверия 0,95 значения критерия анабл гораздо меньше расходится от точного, чем в работах по оценке рисков широко известными ковариационными методами и методами исторического моделирования [2]. При этом точность оценки сохраняется не только в однодневных, но и в более длительных временных горизонтах, что является доказательством высокой адекватности разработанной модели и возможности применения ее для оценивания рисков на финансовых рынках.
Литература
1. БУХБИНДЕР Г. Л., ЧИСТИЛИН К. М. Эмпирическая модель стохастической волатильности финансовых флуктуаций / VI Всероссийская конференция молодых ученых по математическому моделированию: Сборник тезисов. Кемерово, 2005.С. 64-65.
2. ЛОБАНОВ А., ПОРОХ А. Анализ применимости различных моделей расчета Value-at-Risk на российском рынке акций // Рынок ценных бумаг. - 2001. - №2. - С. 65-70.
3. ШИРЯЕВ А. Н. Основы стохастической финансовой математики. М.: Фазис, 2-е изд., 2004 - C. 291-292.
4. CRISAN D., DEL MORAL P., LYONS T. Non-linear filtering using branching and interacting particle systems // Markov Processes Related Fields, Vol. 5, No. 3, 1999. P. 293-319.
5. GORDON N. D., SALMOND SMITH A. F. M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation. // IEEE Proceedings, F-140, 1993. P. 107-113.
6. KARATZAS I., SHREVE S. Brownian motion and stochastic calculus , Springer Verlag, New York, 1991, Second Edition. P. 358-359.
7. PITT M., SHEPHARD N. Filtering via simulation: auxiliary particle filter // Journal of the American Statistical Association, 1999. P. 590-599.
8. STEIN E. M., STEIN J. C. Stock Price Distributions with Stochastic Volatility: An Analytic Approach // Rev. Financial Studies, 1991, 4. P. 727-752.
Статья представлена к публикации членом редакционной коллегии А.И. Орловым.
