Применение методов Монте-Карло для оценки опционов европейского типа

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Методы Монте-Карло можно определить как численные методы решения задач при помощи моделирования случайных величин. Эти методы находят широкое применение в различных областях, в частности, в финансовой математике, например, их можно применить для определения цены опциона европейского типа.

Рассмотрим эволюционирующий в непрерывном времени финансовый рынок, состоящий из двух типов активов: акции (рисковые активы, S) и облигации (безрисковые активы, B). Для моделирования рынка будем использовать модель Блэка–Шоулса. Данная модель после ряда преобразований представляет собой систему стохастических дифференциальных уравнений следующего вида:

$dB_t = rB_tdt, B_0 \gt 0, (1)$

$dS_t = S_t((r-\delta)dt - \sigma dW_t), S_0 \gt 0. (2)$

Здесь первое уравнение представляет собой обыкновенное дифференциальное уравнение, описывающее поведение цены облигации $B_t$, на изменение которой влияет процентная ставка $r$. Уравнение (2) является стохастическим дифференциальным уравнением, описывающим эволюцию цены акции $S_t$. В этом уравнении: $\sigma$ - волатильность акции, $\delta$ - ставка дивиденда, $W_t$ - винеровский случайный процесс.

Сделаем несколько предположений касательно рынка:

1. пусть модель будет без дивидендов: $\delta = 0$;

2. процентные ставки постоянны;

3. если время измеряется в годах, то все рассмотренные процентные ставки участвуют в уравнениях как числа от 0 до 1 и выражают проценты годовых, измеренные в долях.

При этих предположениях рассмотренная система имеет аналитическое решение:

$S_t = S_0 e^{(r - \sigma^2/2)t + \sigma W_t}$.

Опцион, как производный финансовый инструмент, представляет собой контракт между сторонами $P_1$ и $P_2$, который дает право стороне $P_2$ в некоторый момент времени $t$ в будущем купить у стороны $P_1$ или продать стороне $P_1$ акции по зафиксированной в контракте цене $K$. За это право сторона $P_2$ выплачивает фиксированную сумму (премию) $C$ стороне $P_1$. При этом K называется ценой исполнения опциона (страйк, strike price), а $C$ – ценой опциона.

Основная идея заключения контракта состоит в игре двух лиц – $P_1$ и $P_2$. Вторая сторона выплачивает некоторую сумму $C$ и в некоторый момент времени $T$ (срок выплаты, зафиксированный в контракте) принимает решение: покупать акции по цене $K$ у первой стороны или нет. Решение принимается в зависимости от соотношения цены $S_T$ и $K$. Если $S_T \lt K$, покупать акции не выгодно, сторона $P_1$ получила прибыль $C$, а $P_2$ – убыток $C$. В случае, если $S_T \gt K$, сторона $P_2$ покупает у первой акции по цене $K$, в ряде случаев получая прибыль.

Справедливая цена опциона – цена, при которой наблюдается баланс выигрыша/проигрыша каждой из сторон. Определяется она как средний выигрыш стороны $P_2$:

$C = \mathop{\mathbb{E}}(e^{-rT} max\{S_T - K, 0\})$.

Точное решение этой системы уравнений известно:

$C = S_0N(d_1) - Ke^{-rT}N(d_2),$

$d_1 = ln{\frac{S_0}{K}} + \frac{(r + \sigma^2 / 2)T}{\sigma\sqrt{T}},$

$d_2 = d_1 - \sigma\sqrt{T}.$

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

1.2 Математическое описание алгоритма

Решаемая задача имеет следующий вид:

$C = e^{-rT}\int^{+\infty}_{-\infty}{f(z)\varphi(z)dz}, (3)$

$f(z) = max\{S_0 e^{(r - \sigma^2/2)T + \sigma\sqrt{T}z} - K, 0\}, (4)$

$\varphi(z) = \frac{1}{\sqrt{2\pi}} e^{-z^2/2}.$

Для вычисления этого интеграла будем использовать метод Монте-Карло. Если мы рассмотрим случайную величину $X = f(z)$, то математическое ожидание этой случайной величины будет совпадать с формулой (3) с точностью до коэффициента дисконтирования. Формула для вычисления искомого интеграла имеет вид:

$\hat{C} = \frac{e^{-rT}}{N} \sum^{N}_{i=1}{f(z_i)},$

$f(z_i) = max\{S_0 e^{(r - \sigma^2/2)T + \sigma\sqrt{T}z} - K, 0\}, z_i \sim N_{0,1}.$

1.3 Вычислительное ядро алгоритма

Вычислительное ядро алгоритма составляет вычисление следующей суммы:

$\sum^{N}_{i=1}{f(z_i)}.$

Поскольку время, требуемое на операцию сложения в реальных примерах мало по сравнению с временем, затрачиваемом на вычисление функции $f(z_i)$, для производительности вычислений становится маловажным, каким конкретным образом осуществлять сложение и умножение.

1.4 Макроструктура алгоритма

Алгоритм представляет собой сумму значений функции от случайной величины.

1.5 Схема последовательной реализации алгоритма

Формулы для реализации алгоритма приведены выше $(1)$. С помощью генератора случайных чисел в цикле получаются точки, и в них вычисляется соответствующая функция. Результаты складываются, полученная сумма умножается на коэффициент. При последовательной реализации для экономии памяти можно сразу не генерировать массив всех точек, в которых должна быть вычислена функция, а получать их последовательно.

1.6 Последовательная сложность алгоритма

Время, затраченное на вычисление интеграла растет как $Nt$, где $N$ – количество случайных точек, $t$ – время, необходимое для вычисления функции.

1.7 Информационный граф алгоритма

Информационный граф алгоритма представлен на рисунке 1,

Рисунок 1 – Информационный граф алгоритма

где вершина $f$ – вычисление функции $f(z)$, $rnd$ – генерация случайного числа. Анализ этого графа производится в последующих разделах описания алгоритма.

1.8 Ресурс параллелизма алгоритма

Для вычисления интеграла $(3)$ методом Монте-Карло требуется последовательно выполнить следующие ярусы:

1) вычисление случайного числа $z$,

2) вычисление функции $f$ на основании полученного случайного числа,

3) сложение полученных результатов.

В общем случае, все вычисления функции $f$ могут производиться на разных узлах, следовательно, алгоритм идеально поддается распараллеливанию. При классификации по высоте ярусно-параллельной формы этот алгоритм имеет сложность 3. При классификации по ширине – $O(N)$.

1.9 Входные и выходные данные алгоритма

Входные данные для вычисления алгоритма: число точек $N$, волатильность акции $\sigma$, время жизни опциона $T$, страйк $K$, процентная ставка $r$, стоимость акции в начальный момент времени $S_0$. Выходные данные: справедливая цена опциона $C$.

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложностей в случае неограниченных ресурсов является линейным.

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

2 Программная реализация алгоритма

2.1 Особенности последовательной реализации алгоритма

Код последовательной реализации на языке C может выглядеть так:

    double price(int N, double sigma, double T, double K, double S0, double r) {
	coeff = exp(-r*T) / N, result = .0, sigma_coeff = (r - sigma * sigma / 2) * T, sigma_T_coeff = sigma * sqrt(T);
	for (int i = 0; i < N; i++) {
	    result += max(0, S0 * exp(sigma_coeff + sigma_T_coeff * normalrand(i)) - K);
	}
	return coeff * result;
    }

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

2.2 Локальность данных и вычислений

Данные в представленном алгоритме обладают хорошей локальностью, так как на каждом шаге алгоритма приходится использовать значения переменных sigma_coeff и sigma_T_coeff, что повышает вероятность попадания этих значений в кэш.

2.3 Возможные способы и особенности реализации параллельного алгоритма

Один из возможных вариантов реализации алгоритма – независимое вычисление значений функции $f(z)$ на различных вычислительных узлах, например, с помощью технологии MPI. «Подводным камнем» такой реализации может стать генерация случайных чисел, которая должна быть независимой на каждом вычислительном узле. Желательно использовать генератор случайных чисел с достаточно большим периодом.

2.4 Масштабируемость алгоритма и его реализации

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

2.6 Выводы для классов архитектур

При достаточно большом числе точек N возможна реализация алгоритма на GPU для достижения большей скорости выполнения, чем на обычных процессорах.

2.7 Существующие реализации

Есть открытая реализация методов Монте-Карло, являющаяся частью GSL (GNU Scientific Library). Данная библиотека включает в себя как классический алгоритм, так и два адаптивных: «Vegas» и «Miser».

Реализация алгоритма доступна тут: [1], в многокарточном варианте - тут: [2]