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

УДК 519.245 Вестник СПбГУ. Сер. 1. 2013. Вып. 1

ПАРАЛЛЕЛЬНЫЙ МОНТЕ-КАРЛО

МЕТОД ОЦЕНКИ АМЕРИКАНСКИХ ОПЦИОНОВ*

А. В. Дмитриев1, С. М. Ермаков2, К. С. Ермаков3

1. С.-Петербургский государственный университет, аспирант, [email protected]

2. С.-Петербургский государственный университет,

д-р физ.-мат. наук, профессор, [email protected]

3. Независимый консультант в области финансов

и информационных технологий, [email protected]

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

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

Стоимость акций, на которых построен опцион, обозначим через Б(£), а период существования опциона через [0, Т]. Ф. Блэк и М. Шоулз показали, что цена опциона удовлетворяет следующему уравнению в частных производных:

дР а2Б2 д2Р дР -т+ — дБ1 + г8дБ-гР = 0>

где £ — время, Б — цена акции, Р = Р(Б, £) — цена опциона, а — постоянная волатиль-ность, г — безрисковая ставка. Предположения и ограничения, при которых это уравнение имеет место, можно найти в доступной литературе (например, [3]).

Наиболее известными двумя подходами, позволяющими устранить подвижную границу, являются специальная замена переменных и метод штрафных функций [4]. После замены переменных полученное нелинейное уравнение решают методом последовательной линеаризации.

Мы рассмотрим метод штрафных функций совместно с дискретизацией задачи, что дает более простые алгоритмы сравнительно с методом Ньютона.

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00769-а). © А.В.Дмитриев, С.М.Ермаков, X. С. Ермаков, 2013

2. Метод штрафных функций. Суть метода заключается в добавлении к уравнению в частных производных небольшого слагаемого, которое представляет собой непрерывную функцию f (P), нелинейно зависящую от цены опциона P. При этом получается нелинейное дифференциальное уравнение с фиксированной областью определения. Как показано в работе [5], решение этого уравнения будет обладать свойствам аналитического решения задачи нахождения цены американского опциона при достаточной малости добавляемого слагаемого. Рассмотрим уравнение Блэка—Шоулза с такой добавкой,

дР а^^д^Р дРе n m

с граничными условиями

Pe(S,T) = max(K - S, 0), lim Pe(S,t)=0, Pe(0,t) = K, (2)

где

fc

№ = PE + £ — К + S' ^

C > rK — положительная константа. Здесь индексация призвана подчеркнуть зависимость решения уравнения (1)—(2) от f. Справедлива следующая теорема [5].

Теорема 2.1. Пусть P —единственное решение уравнения Блэка—Шоулза, а PF единственное решение уравнения (1) для е > 0. Тогда Ре —> Р в -Loo(Qt) пРи £ ^ 0; где qt = [0, T] х Q, Q С Rm открытое и ограниченное.

Далее индекс f будет опускаться. Пусть 0 < f ^ 1 и dP a2S2 д2P dP

^ + — 9Si+rS--rP + nP)= 0, S<>0, 0<t<T. (4)

Как уже отмечалось, при стремлении f к нулю решение задачи (4), (3), (2) будет стремиться к цене американского опциона. Получившаяся задача является нелинейной с фиксированной областью определения.

Для численного решения рассмотренной задачи (4), (3), (2) могут быть использованы различные методы, их обзор можно найти в [3]. Мы же рассмотрим использование метода конечных разностей.

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

Для того чтобы численно решить систему (4), (3), (2), введем —большое значение S, в котором предполагается выполнение граничного условия на бесконечности, а именно, P(Stt,t) = 0. Для решения задачи (4), (3), (2) методом конечных разностей в области S > 0, 0 < t < T введем сетку

ST М ' N

Si = iAS, для i = 0,..., M, tj = j^t, для j = 0,...,N, Pi,j = P (Si, tj).

Полунеявная разностная схема для уравнения (4) имеет вид

Рг,з+1 - Рг,з , (гАБ)2а Р-^ - 2Р,- + Рт

Аг

+

+ гАБг

Рг+1,3 Ръ—1,3

2АБ

+

= 0, г = 1,...,М - 1, 3 = 0,...^ - 1,

2 (АБ)2

еС

Рг,з+1 +£-К + гАБ

Р, , м = тах(К - гАБ, 0), г = 1,...,М - 1, Рм,з =0, 3 =0,...,Ж - 1, Ро ,3 = к, 3 =0,...,М - 1.

- ГР, , 3 +

После перегруппировки слагаемых получим

г2а2 Аг + ггАг

_ г2cr2Дt-rгДt

^ = 2(1 + г2ст2 Дг + гДг) +

+

2(1 + г2а2Дг + гАг;)Д+и + АгеС

+

(1+ г2а2Аг + гАг) (1+ г2а2Аг + гАг)(Р,,3+1 + е - к + гАБ)' г =1,...,М - 1, 3 =0,...,м - 1.

Введем обозначение Х3 = (Р1,з, Р2,..., Рм-1,3 )Т и заметим, что в общем эта система является нелинейной, однако на каждом шаге по времени получается линейная система вида

Х3 = АХз + ^ (Хз+1), (5)

А

0 С1 0. .0 0 \

а2 0 С2 . .0 0

0 аз 0. .0 0

0 0 0 . .0 см - 2

0 0 0 . ам -1 0 /

(6)

где

г2а2Аг — ггАг

г2а2 Аг + ггАг

г 2(1 + г2а2 Аг + гАг)' г 2(1 + г2а2 Аг + гАг)

а вектор-функция ^ имеет вид

^ (Х3+1) = (/1(Р1,з-+1),/2(Р2,з-+1 ),...,1м-1(Рм-1,3 + 1 ))Т

Р,з+1 АгеС

б, = гАБ, г = 1,...,м- 1,

1г(Р,3+1)

+

(1 + г2а2Аг + гАг) (Р,,3+1 + е - к + гАБ)(1 + г2а2Аг + гАг)' г = 1,...,м-1.

При каждом фиксированном 3 система (5) является системой линейных алгебраических уравнения (Х3+1 вычислено на предыдущем шаге, Х3 —неизвестно). Таким образом, последовательно решив N систем уравнений, начиная с 3, равного N - 1, и до 0, мы найдем цену опциона в начальный момент времени.

3. Стохастические методы. Теперь мы перейдем к основному содержанию работы — применению стохастических методов и исследованию их устойчивости. Докажем следующую лемму.

Лемма 3.1. При каждом фиксированном ] существует Дг > 0 такое, что спектральный радиус .матрицы |А| (6) меньше единицы.

Доказательство. Рассмотрим сумму элементов строки матрицы |А|, состоящей из модулей элементов матрицы А:

, , , , , -гг)\ + г2а2 +гг)

Ы + Ы = г = 1, .. ., М — 1.

При фиксированном ДБ, уменьшая Дг, можно сделать сумму меньше единицы для любого г, сделав норму (в качестве нормы берется максимум из сумм модулей элементов строк) матрицы |А|, а тем самым и ее спектральный радиус, меньше единицы. Можно получить оценку для Дг:

Д£(|(г2ст2 — гъ)\ + г2ст2 + г г) тах(г2ст2, г г)

2(1 + ;2а2Д£ + гД0 = ТТ^дГГТд* < ' ( '

Заметим, что неравенство

г2а2Дг

< 1

1 + г2а2Дг + гДг выполнено для любых Дг и ДБ. Рассмотрим второе неравенство:

тгДг

1 + г2а2 Дг + тДг Напомним, что 1 < г < М, поэтому

< 1.

т%Дг тМДг

< . . о . ■ . , . < 1.

1 + г2а2Дг + тДг - 1 + а2 Дг + тДг После преобразований получим

Дг(т(М —1) - а2) < 1.

При т(М — 1) — а2 < 0 предыдущее неравенство, а следовательно и (7), выполнено. Если же т(М — 1) — а2 > 0, то для выполнения неравенства (7) достаточно следующего условия на Д;:

1

< г(М — I) — а2'

□

Доказанная лемма означает, что можно находить решение системы методом Монте-Карло (возможно также применение квази Монте-Карло). Оба эти метода обладают свойствами естественного параллелизма (ПР-алгоритмы). В этом случае вычисления будут производиться по формуле

& = (Г—А (Мз+1)), (8)

где £3 —оценка Х3, волна в (8) означает, что обращение матрицы и умножение ее на вектор осуществляется посредством метода Монте-Карло, а £3+1 (N3+1) = £3+1 —среднее оценок Х3+1, посчитанных по N3+1 независимым траекториям. При переходе со слоя на слой мы будем иметь дело с двумя видами ошибок: случайной и детерминированной. При решении системы (5) используются несмещенные оценки, однако в силу нелинейности функции Г(Х) будет возникать смещение при переходе со слоя на слой. Возьмем математическое ожидание от правой части (8):

Е ((/-А)-1Г ^3+1(^+1))) =

= Е((Т-А)-1)Е(Г(£3+1 (N3+1))) = (I - А)-1Е(Г(£3+1 (N3+1))).

Первое равенство выполнено в силу того, что на разных слоях используются независимые траектории. Рассмотрим разложение Г(Х) в ряд Тейлора:

Г (£3+1 (N3+1)) =

= р (£ ) = + (^ Е - 1 +

N, + 1

1=1 I \ г=1

I 1 42

Возьмем математическое ожидание от левой и правой части получившегося равенства, предварительно отметив, что Е£3+1 = Е£3+1:

2^3+1

где Г"(Х3+1) —диагональная матрица (М - 1) х (М - 1) с элементами

АгеС

сп __• _ -| Л/Г л

и ~ (1 + г2а2Аг + гДг)(Д,з+1 + е - К + гАБ)3 ' 1 "'' '

В статье [5] было показано, что при Аг < е/С решение уравнения (5) удовлетворяет неравенству

> тах(К - Б\, 0) У г, 3.

Оценим /!':

АгеС Агс

< - < - <

и - (1 + 12а2АЬ + гДг)(тах(^ - 0) + е - К + ¡АБ)3 ~ (1 + 12а2АЬ + гАЬ)е2 ~ 1

~ (1 + ¡2а2 АЬ + гДг)е '

Из этого следует, что при достаточно больших N3 (порядка 1/е и больше) смещением можно будет пренебречь.

4. Условия стохастической устойчивости. Обратимся к рассмотрению случайной ошибки. При переходе по времени со слоя на слой эта ошибка может оставаться ограниченной, а может неограниченно расти с ростом N. Поэтому необходимо исследовать ее поведение. Обозначим £ = ^ — Xз —вектор-столбец случайных ошибок на каждом слое. Выразим еои£з через еои£з+1 в удобной для дальнейшего анализа форме. Из (8) следует

Хз + £э = (1——Л )-\г (Х+ + )). (9)

Напишем разложение функции Г в ряд Тейлора:

+ г '"Х> з

Г X + Ез) = Г X) + Г '(Хз £ + Г ''(Хз )£? + ...

Отметим, что при достаточно малых е и АЬ слагаемым Г''(Xз )£2 и всеми последующими можно будет пренебречь. То есть далее считаем, что Г(Хз + £3) «

Г(Хз) + Г'(Хз)£з. Положим также (I — А) = (I — А)-1 + 6, где Е6 = 0. Равенство (9) перепишется:

Хз + £з = ((I — А)-1 + 6)(Г (Хз+1) + Г '(Хз+1 )£+).

Учитывая, что Хз удовлетворяет (5), получим для £з выражение

£з = 6Г (Хз+1) + (1 — А)-1 Г '(Хз+1 )£з+1 + 6Г '(Хз+1 )£з+1

и соответственно для ¿т, заметив прежде, что Г' —диагональная матрица и, следовательно, (Г')т = Г':

£Т = Г (Хз+1)т 6т + £Т+1 Г '(Хз+1 )((1 — А)-1 )т + ¿^Г '(Х3+1)6т.

Рассмотрим Е(£з,¿т) и введем мультииндекс Ь = (11,12), принимающий (М — 1)2 значений от (1,1) до (М — 1,М — 1). Методами, описанными в [1], можно прийти к

\ть у

(М-1,М-1) ^ м 1

Ь \\ь=(1,1)

^ а- • а- • уз+1

+

-0,-3 = 1

/ у гз ,г2 ± Ь

Ь = (-1 ,-2 )

М-1

Е Е(6го А 6-3,-2 / /4 П+1 • • А ?з+1, (10)

Ь=(-1 ,-2 )

+

го ,гз = 1

где Уэь —элементы матрицы еои£з = Е(£з¿Г), аь,к —элементы матрицы (I — А)-1 Г '(Хз+1), 6го,г16г3гг2 —элементы матрицы 6, /'1 ,/'2 —диагональные элементы матрицы Г'(Хз+1). Матрица Тз+1 объединяет слагаемые, не зависящие от ¿з+1. Теперь вытянем матрицу еои£з в столбец и будем рассматривать далее \\УЬ II как вектор длинны (М — 1)2 , тогда равенство (10) перепишется в форме

\\У1 у = (А + V)\\Yl+1\\ + Тз+1, (11)

А = \\аго,г1 а.гз,г2\\, V = \/'1 Е(6гоЛ6-зг)/2\\ —матрицы (М — 1)2 х (М — 1)2. При этом собственные векторы А суть угц>т, г,] = 1,...,М — 1, где уг — собственный вектор-

столбец матрицы (I — А) 1Г', а собственные числа А суть ХгХз, г,] = 1,...,М —

1, где Хг — собственное число матрицы (I - А) 1Г'. Матрица Р представляет собой произведение МГ, где Г —диагональная матрица с диагональю ||/'/' У^НС?^) а элементы матрицы М суть ковариации вектора, составленного из элементов матрицы погрешностей 6. Элементы матрицы М имеют порядок 1/Nз, где N3 —количество независимых траекторий на слое с номером 3, то есть М = 1/NзМ'. Оценим /', выражение для которого имеет вид

*> =_1_Л___

и 1 + ^а2Аг + гАг\ (Рг,з+1 + £ - к + ¡АБ)2

найдем условие на Аг, при котором разность, стоящая в скобках, больше нуля,

АгеС _ АгеС

1

(Рг,3+1 + е - К + гАБ)2 " откуда получим условие

(е - К)2

1

(е - К)

2>

0,

Аг

еС

При Аг из (12) норма матрицы Г' меньше единицы. Таким образом,

\\¥[\\=(Л+— МР)\\УЬ

3+1

(12)

(13)

и матрица ковариаций соу£ 3 будет оставаться ограниченной, если норма ||Д + 1/NзМ'Г|| будет меньше единицы и будет расти экспоненциально, если норма ||Д + 1/NзМ-Р|| будет больше единицы. Вышесказанное и выражение (13) влекут следующие утверждения.

Утверждение 4.1. Для стохастической устойчивости алгоритма (8) необходимо и достаточно, чтобы максимум из модулей первых собственных чисел матриц Л+ -^-ЛА'Р был меньше единицы.

Утверждение 4.2. Если первые собственные числа матриц (I - А)- 1Г' по модулю строго меньше единицы, то существует такое натуральное число N0, что при всех N 3 > N0 алгоритм (8) будет стохастически устойчивым.

Необходимо оценить первое собственное число матрицы (I-А)- 1Г'+1, где Г'+1 = Г'(Х3+1). Матрицу Г3+1 можно представить в виде произведения диагональных матриц Г1Г2 с элементами, равными соответственно 1 АгеС

1 + г2а2Аг + гАг

и1

(Рг,3+1 + е - К + гАБ)2

1,...,М - 1.

Далее,

где

(I - А)-1Г3+1 = (Г1 (I - А)) Г = (I - АгА)-1 (I - АгГ)

а4:

/Ъ1 с1 0. .0 0 \

а2 Ъ2 с2 . .0 0

0 аз Ъз . .0 0

0 0 0 . Ъм-1 см - -2

0 0 0 . ам-1 Ъм- 1

1/1 0 0 /2

Г

0

0

/з

/м-2 0 0 /м -1/

г2а2 - т"1 22 г2а2 + гг сч = -^-' Н = -г а - г, с» = ---,

еС

(Дз+1 + е - К + гД6)2 При достаточно малых Дг можно считать (I — ДгЛ)-1 « I+ ДгЛ, следовательно,

(I — ДЛ)-1 (I — Д] « (I + дгЛ)(1 — Д] «I + (Л — Р)Дь.

Если матрица Л — ] отрицательно определенная, то при достаточно малых Дг первое собственное число матрицы I + Дг(Л — ]), а следовательно и матрицы (I — Л)-1 ], меньше единицы. Можно показать, что при г > а2/3 матрица Л — ] отрицательно определенная. Также выполнение условия

~ (4г2 + 4г + 1)<т4 — (12г2 + 12г + 6)<т2г — Зг2 " > 4(г+ 1)2а2 +4г2а2 +4г ( ^

для любого г обеспечит отрицательную определенность. Это условие также легко проверить в процессе вычислений. Показать отрицательную определенность матрицы Л — ] в общем случае не удалось. Более того, в некоторых, далеких от реальных ситуаций случаях возможно нарушение отрицательной определенности. Тем не менее, в большинстве случаев стохастическая устойчивость имеет место и это означает, что решение задачи можно поручить разным процессорам. Окончательный результат получается осреднением результатов отдельных процессоров. Кроме того, можно оценить ошибку, вычисляя дисперсию оценок.

5. Численные эксперименты. Ниже приведены результаты решения систем (5) методом Монте-Карло и детерминированным методом. Расчеты производились для следующих параметров: Т = 0.75, = 100, М = 100, N = 700, а = 0.15, г = 0.055, К = 35, е = 0.001, С = 2. На графиках изображена цена опциона в зависимости от цены акции в момент времени г = 0 (рис. 1).

Количество моделируемых траекторий для оценки каждого узла сетки 50 — для первого графика, 100 —для второго. Сплошная линия соответствует решению, полученному методом Монте-Карло (нанесены средние значения оценок). Доверительный интервал при уровне доверия а = 0.95 достаточно мал и на графиках не изображен. Штриховая линия соответствует решению, точно посчитанному детерминированным методом.

Также были проведены эксперименты для следующих параметров: Т = 0.75, ^ = 100, М = 100, N = 900, а = 0.7, г = 10-4, К = 35, е = 10-4, С = 2 (рис. 2). Количество моделируемых траекторий для оценки каждого узла сетки: 200 — для первого графика, 250 — для второго. Сплошная линия соответствует решению, полученному методом Монте-Карло, штриховая — точно посчитанному детерминированным методом.

Условие (14) очевидным образом может быть наиболее ограничительным при малых г. Случай г = 0 не реален, но интересно было проследить поведение алгоритма при малых г.

Следующие примеры (г = 10-4) обнаруживают нарушение устойчивости (даже детерминированной!) в начальной области (наличие «горба» на графиках) и сравнительное увеличение статистической ошибки.

Рис. 1. Зависимость цены опциона от цены акции для разного количества моделируемых траекторий: 50 — слева; 100 — справа.

Рис. 2. Зависимость цены опциона от цены акции для разного количества моделируемых траекторий: 200 — слева; 250 — справа.

Заключение. Проведенное исследование показывает, что при выполнении условий устойчивости (они, как правило, выполняются при умеренных значениях N) возможно применение метода Монте-Карло. Вычисления могут проводиться на многих независимых процессорах (облачные вычисления). Мы рассмотрели сравнительно простую модель. Можно ожидать, что полученные выводы будут справедливы для более сложных моделей рынка, когда применение метода Монте-Карло является особенно актуальным.

Литература

1. Ермаков С. М. Метод Монте-Карло в вычислительной математике (Вводный курс). СПб.: Бином. Лаборатория знаний, 2009. 192 с.

2. Ермаков С. М. Параметрически разделимые алгоритмы // Вестн. С.-Петерб. ун-та. Сер. 1. Вып. 4. 2010. С. 68-73.

3. Duffy D. Finite difference methods in financial engineering: a partial differential equation approach. John Wiley & Sons Ltd, 2006. 423 p.

4. Zvan R., Forsyth P. A., Vetzal K. R. Penalty methods for American options with stochastic volatility // Journal of Computational and Applied Mathematics. no. 218. P. 91-199.

5. Nielson B. F., Skavhaug O., Tvelto A. Penalty and front-fixing methods for the numerical solution of American option problems //J. Comp. Finan. no. 4, 2002.

Статья поступила в редакцию 20 сентября 2012 г.