Санкт-Петербургский государственный университет
Кафедра компьютерного моделирования и
многопроцессорных систем
Войцеховская Виктория Анатольевна
Выпускная квалификационная работа бакалавра
Оценка стоимости европейского опциона
методом Монте-Карло
Направление 010300
Фундаментальная информатика и информационные технологии
Научный руководитель,
д-р физ.-мат. наук,
профессор
Богданов А.В.
Санкт-Петербург
2016
Содержание
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Предметная область . . . . . . . . . . . . . . . . . . . . . . . . . .
6
Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
. . . . . . . . . . . . . . . . . . . . . . . . . . .
9
Последовательные алгоритмы оценки опционов . . . . . .
10
Введение
Обзор литературы
1
2
1.1
1.2
. . . . . . . . . .
Вычислительные мощности . . . . . . . . . . . . . . . . . . .
Параллельные вычисления в MATLAB . . . . . . . . . . . .
2.2.1 О MATLAB и Parallel Computing Toolbox . . . . . . .
2.2.2 Параллельный алгоритм в MATLAB . . . . . . . . .
Технология параллельных вычислений OpenCL . . . . . . .
2.3.1 Описание технологии OpenCL . . . . . . . . . . . . .
2.3.2 Генератор случайных чисел . . . . . . . . . . . . . . .
2.3.3 Редукция выходного массива . . . . . . . . . . . . . .
2.3.4 Структура программы на OpenCL . . . . . . . . . . .
2.3.5 Параллельный алгоритм расчета опционов на OpenCL
Параллелизм для ускорения вычислений
2.1
2.2
2.3
3
Численные методы оценки опционов . . . . . . . . . . . . . .
Алгоритмы оценки опционов, основанные на методе МонтеКарло . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Стохастическое дифференциальное уравнение . . . .
1.2.2 Континуальный интеграл . . . . . . . . . . . . . . . .
1.2.3 Дифференциальное уравнение в частных производных
Анализ полученного ускорения вычислений . . . . . . . . .
3.1
3.2
3.3
Ускорение в MATLAB . . . . . . . . . . . . . . . . . . . . . .
Ускорение при помощи OpenCL . . . . . . . . . . . . . . . .
Сравнение времени работы реализаций в MATLAB, на языке
C и OpenCL . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
10
11
16
20
24
24
25
25
26
27
27
29
31
33
35
37
37
40
46
Выводы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
. . . . . . . . . . . . . . . . . . . . . . . . . .
51
Список литературы
2
Введение
Метод Монте-Карло - метод численного решения математических задач
с помощью моделирования большого количества случайных величин и нахождения их математического ожидания. Вычисление интегралов является
основной направленностью этого метода, так как математическое ожидание непрерывной случайной величины выражается через интеграл.[1][2]
Метод Монте-Карло имеет применение во многих областях, как то:
∙ задачи финансовой математики - моделирование рыночных ситуаций;
∙ задачи теории игр;
∙ производственные задачи - моделирование сложных систем и сетей;
∙ задачи ядерной и статической физики, и другие.
Популярность этому методу принесла его существенная простота. Данный метод основан на теории вероятности, а конкретнее, на Центральной
предельной теореме. Таким образом, ошибка вычислений, проводимых на
основе этого метода, сильно зависит от количества моделируемых испытаний. Чем больше испытаний мы проводим, тем меньше ошибка, и в то
же время больше временные затраты на расчеты. Тем не менее, для некоторых задач, которые требуют высокой точности вычислений, даже достаточно большое количество испытаний, проведенных в рамках метода
Монте-Карло, не даст желательной точности.
Однако мы будем рассматривать применение этого метода в рамках
моделирования экономических ситуаций. Обычно в таких задачах требуется спрогнозировать примерную траекторию развития рыночной ситуации. Например, нам не нужно знать сколько будет стоить акция через год
с точностью до копейки, но нам нужно знать, будет ли цена на акцию
расти или падать, и чем быстрее мы это узнаем, тем лучше. Таким образом, для прогнозирования в рамках экономической и финансовой области
в приоритете оказывается скорость вычислений, нежели их точность. И
тут оказывается полезным свойство хорошей распараллеливаемости метода Монте-Карло, путем распределения численных статических испытаний
3
по отдельным процессорам.
В настоящее время большую популярность приобрели различные технологии параллельных вычислений, в том числе с использованием GPGPU
(General Purpose GPU). Это связано с постоянно растущей сложностью
актуальных задач, а так же с растущим объемом данных для обработки. Суть таких вычислений в том, чтобы части программы выполнялись
одновременно в наборе потоков, взаимодействующих друг с другом. Использование параллельных компьютеров (компьютерных систем с набором
процессоров, работающих одновременно) обосновано идеей о том, что если
одному процессору на выполнение задачи требуется 𝑝 времени, то 𝑛 процессорам на это должно потребоваться в 𝑛 раз времени меньше. Но, конечно,
это идеальный случай, и в реальности такого добиться не удастся. Однако,
ускорить работу вполне возможно, в той или иной степени, в зависимости
от алгоритма и имеющейся аппаратной системы.[3]
Примерами технологий для параллельных вычислений могут быть
OpenMP, MPI, CUDA, OpenCL. Так же в пакете прикладных программ
MATLAB существует набор средств для написания параллельных алгоритмов, объединенных в MATLAB Parallel Computing Toolbox.
В данной работе будут исследоваться алгоритмы оценки европейского
опциона, как примера классического вида, и азиатского, в качестве усложнения первоначальной задачи. Рассматриваемые алгоритмы, основанные
на методе Монте-Карло, будут реализованы последовательно и параллельно с помощью технологий MATLAB Parallel Computing Toolbox и технологии гетерогенных вычислений OpenCL и языка C. Анализ скорости работы
разных подходов и технологий будет проведен на высокопроизводительном
кластере Ресурсного Центра "Вычислительный Центр СПбГУ"[4].
Опцион - вид экономических деривативов, который является достаточно гибким инструментом работы на бирже и позволяет снизить риски манипулирования активами. Подробнее опционы будут описаны в главе "Предметная область". Большую роль в теории оценки опционов сыграла статья
4
"The pricing of options and corporate liabilities" 1973 года (Myron Scholes,
Fischer Black)[5], в которой была описана модель Блэка-Шоулза, давшая
аналитическую формулу для классических, или ванильных, опционов. Но
для опционов другого вида, экзотических, формулы нет, поэтому приходится использовать численные методы для приближенного нахождения цены
таких опционов, что и приводит нас к методу Монте-Карло.
Таким образом, в данной работе будет исследован вопрос об оценке опционов методом Монте-Карло и ускорении соответствующих вычислений
с помощью программных средств MATLAB Parallel Computing Toolbox и
технологии гибридных вычислений OpenCL.
5
Предметная область
Опционы
Опцион - один из видов деривативов; контракт между покупателем и
продавцом, который дает потенциальному покупателю или продавцу право,
но не обязательство, купить базовый актив (которым может быть, например, некоторая акция) по заранее обговоренной цене в некоторый момент
в будущем.[6]
Соответственно, выделяют два вида опционов: Call и Put. Call-опцион
дает право покупателю данного дериватива купить базовый актив, а Putопцион - продать.
Цена продажи базового актива, указанная в договоре, называется Strikeценой. Дата, когда должен быть выполнен договор, называется днем истечения.
Существует много разных типов опционов, самые известные из них это
европейский и американский, которые называются стандартными или ванильными опционами. Они заранее оговаривают все детали договора. Но
развитие рынка привело к появлению опционов других видов, так называемых, "экзотических" опционов. Среди них находятся такие опционы, как:
азиатский, барьерный, бинарный, свопцион и др.
В данной работе будут рассмотрены европейский и азиатский опционы.
Европейский опцион является классическим, в нем оговорена Strike-цена и
срок жизни опциона. Такой опцион исполняется в день истечения договора
по оговоренной цене. Стоимость такого Call-опциона можно вычислить по
формуле:
Φ(𝑆) = (𝑆(𝑇 ) − 𝐾)+ ,
(1)
где S(t) - цена акции в момент t, T - дата истечения договора, K - Strikeцена опциона.
Азиатский опцион исполняется по средневзвешенной цене акции за весь
6
период жизни опциона с момента покупки, поэтому формула для данного
случая будет отличаться от формулы (1) тем, что вместо 𝑆(𝑇 ) будет 𝑆𝑎𝑣𝑔 (𝑇 )
- среднее арифметическое цен на актив за период времени (0; 𝑇 ).
Таким образом, вычисление стоимости опциона сводится к прогнозированию цены на базовый актив.
Модель ценообразования опционов Блэка-Шоулза
Модель Блэка-Шоулза[5] - одна из самых популярных моделей оценки
стоимости опционов. Подразумевая, что если опцион торгуется на рынке,
то сам рынок неявным образом устанавливает на него цену, эта модель дает
аналитическую формулу для вычисления стоимости европейского опциона.
Конечно, как и любая математическая модель, она устанавливает некоторые допущения:
∙ дивиденды по базовому активу не выплачиваются во время всего срока действия опциона;
∙ нет дополнительных затрат, связанных с покупкой или продажей актива или опциона;
∙ безрисковая процентная ставка известна и является постоянной во
время всего срока действия опциона;
∙ продажа базового актива ведется непрерывно во время действия опциона, а его цена изменяется согласно модели геометрического броуновского движения;
∙ и др.
Таким образом, модель Блэка-Шоулза часто используется на практике
и является удобной аналитической формулой, но не может оценить более
сложных опционов, чем европейский.
7
Постановка задачи
В рамках данной работы необходимо выполнить следующие задачи:
∙ Исследовать различные подходы к оценке стоимости европейского
(в качестве представителя классических или ванильных опционов)
и азиатского (в качестве представителя экзотических опционов) опциона, основанные на методе Монте-Карло и модели Блэка-Шоулза;
∙ Составить алгоритмы, реализующие исследуемые подходы;
∙ Проверить работу написанных алгоритмов в среде MATLAB;
∙ Написать параллельные алгоритмы для ускорения вычислений в среде MATLAB, проанализировать скорость работы алгоритмов при разных количествах потоков;
∙ Написать параллельные алгоритмы для ускорения вычислений для
технологии OpenCL, проанализировать скорость работы;
∙ Сделать выводы о полученных результатах.
8
Обзор литературы
Источник [6] является замечательным пособием, раскрывающим тему
опционов, их ценообразования, используемых в данной области методов, а
так же дает развернутое представление о расчете стоимости европейского
и азиатского опциона методом Монте-Карло из стохастического дифференциального уравнения.
Статья в двух частях [7] и [8] подробно расписывает применение метода
Монте-Карло к оценке европейского и азиатского опциона на основе стохастического дифференциального уравнения и приводит пример расчетов
в Excel.
Работы [9] и [10] дают понятие о численном интегрировании методом
Монте-Карло, что используется во втором рассматриваемом нами подходе
для нахождения цены опционов различных видов.
Работа [11] иллюстрирует подход к оценке стоимости опционов методом
Монте-Карло, основанном на уравнении в частных производных.
В работе [12] рассматриваются опционы в целом и различные подходы
к оценке их стоимости, в том числе не только метод Блэка-Шоулза, но и
другие.
Для работы с MATLAB очень полезен сайт [13], на котором можно
найти множество авторских руководств, описаний функций и прикладных
пакетов.
Для изучения OpenCL были использованы источники [14] и [15]. В первом из них подробно и с примерами описаны нюансы работы с OpenCL
и организация параллельного приложения. Во втором подробнее рассматриваются подходы программирования на OpenCL, например, MapReduce,
который будет использован в этой работе.
9
1
Последовательные алгоритмы оценки опционов
1.1
Численные методы оценки опционов
Рассмотрим три вида численных методов, используемых для оценки
опционов[6]:
1. Модели конечных разностей - используются, когда динамику цены опциона можно описать в терминах дифференциального уравнения в частных производных.
2. Биномиальные методы - берут за основу модель Блэка-Шоулза,
но являются более гибкими, и поэтому более точными. Неплохо вычисляют цену как европейского, так и американского опциона. По
сути являются частным случаем явной модели конечных разностей.
Но не вычисляет цену экзотических опционов.
3. Методы Монте-Карло - основаны на множественном генерировании случайных величин и вычислении на их основе цены актива в
нужный момент времени. Справляется как с классическими, так и с
экзотическими опционами, является достаточно простым, но требует
много вычислений.
Таким образом, для решения задачи оценки стоимости экзотических, и
в частном случае, азиатских, опционов, лучшим решением будет использовать метод Монте-Карло. Далее мы рассмотрим алгоритмы, основанные
на этом методе, и ускорим их работу.
1.2
Алгоритмы оценки опционов, основанные на методе Монте-Карло
Существует несколько подходов к решению задачи ценообразования
опциона, использующих метод Монте-Карло.[7][8][6] Рассмотрим методы,
основанные на стохастическом дифференциальном уравнении динамики
цены базового актива(далее иногда - СДУ), на континуальном интеграле(далее иногда - КИ) и на уравнении в частных производных(далее иногда
10
- УвЧП) для европейского и азиатского опциона. Оценка стоимости европейского опциона покажет погрешность методов, так как есть возможность
сравнить полученные значения со стоимостью по аналитической формуле
модели Блэка-Шоулза. [5]
Данная формула выглядит следующим образом для Call-опциона:
𝐶(𝑆0 , 0) = 𝑆0 𝑁 (𝑑1 ) − 𝐾 exp−𝑟𝑇 𝑁 (𝑑2 ),
𝑆0
) + (𝑟 + 𝜎 2 /2)𝑇
𝐾
√
,
𝑑1 =
𝜎 𝑇
𝑆0
ln( ) + (𝑟 − 𝜎 2 /2)𝑇
𝐾
√
𝑑2 =
,
𝜎 𝑇
где 𝑆0 - цена базового актива в начальный момент времени, 𝐾 - Strike-цена,
𝑟 - безрисковая ставка, 𝜎 - волатильность, 𝑇 - остаток времени до истечения договора (как часть от года). 𝑁 (𝑥) - функция распределения вероятностей стандартного нормального закона, иными словами, рассчитывающая
вероятность попадания в интервал (−∞, 𝑥) величины, распределенной по
стандартному нормальному распределению.
ln(
Для расчетов и сравнения методов возьмем следующие начальные данные: 𝑆0 = 100, 𝐾 = 110, 𝑟 = 0.1, 𝜎 = 0.3, 𝑇 = 1. Стоимость европейского
Call-опциона 𝐶 в данном случае по аналитической формуле Блэка-Шоулза
равна 12.1310. С этим результатом и будем сравнивать работу алгоритмов.
1.2.1
Стохастическое дифференциальное уравнение
Согласно модели Блэка-Шоулза стохастическое дифференциальное
уравнение динамики цены базового актива выглядит следующим образом:
𝑑𝑆(𝑡) = 𝑟𝑆(𝑡)𝑑𝑡 + 𝜎𝑆(𝑡)𝑑𝐵(𝑡),
11
где 𝐵(𝑡) - стандартный винеровский процесс, 𝐵(𝑡) ∼ 𝑁 (0, 1). Тогда для
европейского опциона справедлива формула:
√
𝜎2
𝑆(𝑇 ) = 𝑆0 exp {(𝑟 − )𝑇 + 𝜎 𝑇 𝜖},
2
(2)
где 𝜖 - это случайная величина, распределенная по стандартному нормальному закону. Рассчитав цену базового актива в день исполнения, не вызывает труда вычислить стоимость опциона по функции выигрыша:
𝐶(𝑆0 , 𝑇 ) = 𝑚𝑎𝑥(𝑆(𝑇 ) − 𝐾, 0)
(3)
Для азиатского же опциона требуется отслеживать не только начальную и конечную цену актива, но и цены на всем временном промежутке
действия опциона. Тогда в каждый момент времени цена будет зависеть от
цены на предыдущем шаге, и формула примет вид:
√
𝜎2
𝑆(𝑡𝑖+1 ) = 𝑆(𝑡𝑖 ) exp {(𝑟 − )𝛥𝑡 + 𝜎 𝛥𝑡𝜖},
2
где 𝛥𝑡 - промежуток времени между двумя измерениями цены. Для каждого набора из 𝑚 значений цен на актив в некоторый момент времени найдем
средневзвешенное значение. И тогда по формуле 𝐶(𝑆0 , 𝑇 ) = 𝑚𝑎𝑥(𝑆𝑎𝑣𝑔 (𝑇 )−
𝐾, 0) находим стоимость азиатского опциона. После этого проводим такие
же вычисления еще 𝑛 раз, и дисконтированное среднее значение цены будет
искомой стоимостью.
Европейский опцион.
Формальный алгоритм для вычисления цены ев-
ропейского опциона:
1. Генерация 𝑛 случайных, нормально распределенных чисел 𝜖 и вы-
числение для каждого из них цены актива 𝑆(𝑇 ) по формуле (2).
2. Нахождение для каждой найденной цены актива 𝑆(𝑇 ) цены опциона
𝐶(𝑆0 , 𝑇 ) по формуле (3).
3. Вычисление среднего арифметического всех полученных цен опциона
𝐶𝑚𝑖𝑑 .
4. Вычисление дисконтированного среднего значения цены опциона по
12
формуле:
𝐶 = 𝑒−𝑟𝑇 𝐶𝑚𝑖𝑑
Азиатский опцион.
Так как азиатский опцион выполняется по сред-
невзвешенной цене актива за всю жизнь договора, алгоритм немного видоизменяется. Будем разделять время жизни опциона на столько частей,
сколько раз будет наблюдаться цена актива. Пусть временной промежуток
𝑇
. Формальный алгоритм для вычис𝑚
ления цены азиатского опциона принимает следующий вид:
1. Генерация 𝑚 случайных, нормально распределенных чисел 𝜖 и вычисление для каждого из них цены актива 𝑆(𝑡𝑖 ) по формуле:
разделен на 𝑚 частей, тогда 𝛥𝑡 =
√
𝜎2
𝑆(𝑡𝑖+1 ) = 𝑆(𝑡𝑖 ) exp {(𝑟 − )𝛥𝑡 + 𝜎 𝛥𝑡𝜖}
2
2. Вычисление среднего арифметического 𝑆𝑎𝑣𝑔 (𝑇 ) всех 𝑆(𝑡𝑖 ) и нахож-
дение цены опциона 𝐶(𝑆0 , 𝑇 ) по формуле 𝐶(𝑆0 , 𝑇 ) = 𝑚𝑎𝑥(𝑆𝑎𝑣𝑔 (𝑇 ) − 𝐾, 0)
3. Повторение пунктов 1-2 𝑛 раз.
4. Вычисление среднего арифметического от n полученных значений
цены опциона и затем дисконтированного среднего значения по формуле:
𝐶 = 𝑒−𝑟𝑇 𝐶𝑚𝑖𝑑
Реализация в MATLAB, сходимость.
В процессе рассмотрения по-
ставленной задачи были реализованы описанные выше алгоритмы в среде
MATLAB(версия R2013b).
Рассмотрим сходимость данного метода. В случае европейского опциона (Рис.1, Рис.2, Рис.3), видно, что метод сходится к аналитическому решению, не смотря на то, что количество генерируемых случайных величин
достаточно мало. Для анализа скорости последовательных и параллельных
алгоритмов будет взято гораздо большее количество, от 1 до 50 миллионов.
13
Рис. 1: По горизонтали N, по вертикали стоимость опциона. Синий цвет метод Монте-Карло, красный - аналитический метод.
Рис. 2: По горизонтали N, по вертикали ошибка вычислений по сравнению
со значением, полученным с помощью аналитической формулы.
14
Рис. 3: По горизонтали N, по вертикали стоимость опциона. Синий цвет метод Монте-Карло, красный - аналитический метод.
Для проведения данного анализа было произведено 300 вычислений с шагом в 200(Рис.1 и Рис.2). На Рис.3 изображен график, полученный в результате проведения 1000 вычислений с шагом в 100.
В случае азиатского опциона, нет аналитической формулы, но можно
всё равно проследить сходимость и сравнить, например, со средним значением получаемых цен(Рис.4).
15
Рис. 4: По горизонтали N, по вертикали стоимость опциона. Синий цвет
- метод Монте-Карло, красный - среднее значение полученных в процессе
стоимостей.
1.2.2
Континуальный интеграл
Решение уравнения Блэка-Шоулза в частных производных может быть
найдено следующим образом:
𝐶(𝑆, 𝑇 ) = exp−𝑟𝑇 𝐸[𝐹 (𝑆𝑇 )],
где 𝐸[·] - математическое ожидание, а 𝐹 (𝑆𝑇 ) - функция выигрыша опциона. Как известно, математическое ожидание можно представить в виде
интеграла, в данном случае такого:
∫︁+∞ ∫︁+∞
𝐸[𝐹 (𝑆𝑇 )] =
...
𝑑𝑥1 . . . 𝑑𝑥𝑛 𝑑𝑥𝑛+1 max(exp 𝑥𝑛+1 − 𝐾, 0)×
−∞
−∞
16
𝑛+1
∑︁
1
1
[𝑥𝑛 − (𝑥𝑛−1 + 𝜇𝛥𝑡)]2 ]}
×
exp{
[−
2
2
(𝑛+1)/2
2𝜎 𝛥𝑡
(2𝜋𝜎 𝛥𝑡)
𝑖=1
Тогда можно вычислить его значение с помощью метода Монте-Карло.
Европейский опцион.
Исходя из выше сказанного, формальный алго-
ритм для расчета цены европейского опциона выглядит следующим образом:
1. Генерация 𝑛 случайных чисел 𝑥𝑖 таких, что 𝑥0 = ln 𝑆0 , а плотность
распределения 𝑥𝑖 зависит от 𝑥𝑖−1 следующим образом:
√
1
2𝜋𝜎 2 𝛥𝑡
exp {−
1
[𝑥𝑖 − (𝑥𝑖−1 + 𝜇𝛥𝑡)]2 }
2
2𝜎 𝛥𝑡
2. Вычисление функции 𝐶(𝑇 ) = max(𝑒𝑥𝑛 − 𝐾, 0).
3. Повторение шагов 1-2 𝑚 раз, затем вычисление среднего значения
полученных цен 𝐶𝑚𝑖𝑑 .
4. Вычисление дисконтированного среднего значения цены опциона по
формуле:
𝐶 = 𝑒−𝑟𝑇 𝐶𝑚𝑖𝑑
Азиатский опцион.
Алгоритм же для азиатского опциона может быть
записан следующим образом:
1. Генерация 𝑛 случайных чисел 𝑥𝑖 таких, что 𝑥0 = ln 𝑆0 , а плотность
распределения 𝑥𝑖 зависит от 𝑥𝑖−1 следующим образом:
√
1
2𝜋𝜎 2 𝛥𝑡
exp {−
1
[𝑥𝑖 − (𝑥𝑖−1 + 𝜇𝛥𝑡)]2 }
2
2𝜎 𝛥𝑡
2. Вычисление среднего значения 𝑥𝑚𝑖𝑑 =
𝐾, 0).
Σ𝑥𝑖
и функции 𝐶(𝑇 ) = max(𝑒𝑥𝑚𝑖𝑑 −
𝑛
3. Повторение шагов 1-2 𝑚 раз, затем вычисление среднего значения
полученных цен 𝐶𝑚𝑖𝑑 .
17
4. Вычисление дисконтированного среднего значения цены опциона по
формуле:
𝐶 = 𝑒−𝑟𝑇 𝐶𝑚𝑖𝑑
Реализация в MATLAB, сходимость.
Реализации методов были осу-
ществлены в среде MATLAB. Рассмотрим графики, полученные в результате анализа работы реализации данного алгоритма для европейского опциона (Рис.5, Рис.6).
Как видно по графикам, метод сходится к аналитическому значению,
ошибка с увеличением числа 𝑛 уменьшается. В случае азиатского опциона
имеем графики, которые показаны на Рис.7 и Рис.8. Так же прослеживается, как метод сходится.
Рис. 5: По горизонтали N, по вертикали стоимость опциона. Синий цвет метод Монте-Карло, красный - аналитический метод.
18
Рис. 6: По горизонтали N, по вертикали ошибка вычислений по сравнению
со значением, полученным с помощью аналитической формулы.
Рис. 7: По горизонтали N, по вертикали стоимость опциона. Синий цвет метод Монте-Карло, красный - среднее значение всех вычисленных цен в
течение испытания.
19
Рис. 8: По горизонтали N, по вертикали ошибка вычислений по сравнению
со средним значением.
1.2.3
Дифференциальное уравнение в частных производных
Европейский опцион.
Для европейского опциона метод, основанный на
дифференциальном уравнении в частных производных, вырождается в метод, основанный на континуальном интеграле ввиду одномерности задачи.
Азиатский опцион.
Метод сводится к методу блуждания по сетке и
затем решению краевой задачи. Берем на области сетку, определяем начальную точку, а затем генерируем вероятности перехода в соседние точки.
Таким образом находим путь, который приводит нас к одной из границ. В
зависимости от границы и конкретной конечной точки, вычисляем искомое значение цены опциона. Формальный алгоритм состоит из следующих
шагов:
1. Разбивание прямоугольника [0, 𝑆𝑚𝑎𝑥 ] × [0, 𝐴𝑚𝑎𝑥 ] сеткой на маленькие
прямоугольники.
2. Выбор начальной точки, в которой требуется найти цену опциона. В
данном случае перейти можно только в три точки. Левая от текущей точки
20
(i-1,j), верхнюю - (i, j+1), а правую - (i+1,j). Вероятности перехода в эти
точки назовем соответственно 𝑝𝑖−1,𝑗 , 𝑝𝑖,𝑗+1 𝑝𝑖+1,𝑗 .
3. Построение пути по следующим правилам:
1. Вычисление на каждом шаге вероятностей перехода в соседние точки
по формулам:
Γ𝑖−1,𝑗,𝑛
𝜎 2 𝑆𝑖 2
=
2(𝛥𝑆)2
Γ𝑖,𝑗+1,𝑛 =
Γ𝑖+1,𝑗,𝑛
𝑄𝑖,𝑗,𝑛 =
𝑆𝑖
𝛥𝐴
𝜎 2 𝑆𝑖 2
𝑟𝑆𝑖
=
+
2(𝛥𝑆)2 𝛥𝑆
1
Γ𝑖−1,𝑗,𝑛 + Γ𝑖,𝑗+1,𝑛 + Γ𝑖+1,𝑗,𝑛
= 𝛥𝑡
𝑝( 𝑖 − 1, 𝑗, 𝑛) = Γ𝑖−1,𝑗,𝑛 · 𝑄𝑖,𝑗,𝑛
𝑝( 𝑖, 𝑗 + 1, 𝑛) = Γ𝑖,𝑗+1,𝑛 · 𝑄𝑖,𝑗,𝑛
𝑝( 𝑖 + 1, 𝑗, 𝑛) = Γ𝑖+1,𝑗,𝑛 · 𝑄𝑖,𝑗,𝑛
При этом, если сумма всех 𝑄𝑖,𝑗,𝑛 = 𝛥𝑡, вычисленных во время построения пути, превысила T, то нужно прервать вычисление. В таком
случае значение 𝑉𝑘 вычисляем по формуле:
𝑉 (𝑆, 𝐴, 0) = 𝑒−𝑟𝑇 max
𝐴𝑗
− 𝐾, 0
𝑇
2. Если была достигнута какая-либо из границ области, то вычисляем
значение 𝑉𝑖 с помощью соответствующих краевых условий:
𝑉 (𝑆𝑚𝑎𝑥 , 𝐴, 𝑡) = max(
𝐴𝑗
𝑆𝑚𝑎𝑥
− 𝐾, 0) +
* (𝑒−𝑟*𝑄𝑠 𝑢𝑚 − 1)
𝑇
𝑟*𝑇
𝐴𝑚𝑎𝑥
𝑆𝑖
− 𝐾, 0) +
* (𝑒−𝑟*𝑄𝑠 𝑢𝑚 − 1)
𝑇
𝑟*𝑇
𝐴𝑗
𝑉 (0, 𝐴, 𝑡) = max( − 𝐾, 0)
𝑇
𝑉 (𝑆, 𝐴𝑚𝑎𝑥 , 𝑡) = max(
4. Таким образом строится 𝑁 различных путей и, соответственно, вычис-
ляется 𝑁 различных значений 𝑉𝑖 . Находим среднее арифметическое всех
21
этих значений, и тогда цена опциона определяется формулой:
𝐶 = 𝑒−𝑟𝑇 𝑉𝑚𝑖𝑑
Реализация в MATLAB, сходимость.
Алгоритм был реализован в
среде MATLAB.
Графики, получившиеся в результате анализа работы программы представлены ниже (Рис.9, Рис.10):
Рис. 9: По горизонтали N, по вертикали стоимость опциона. Синий цвет метод Монте-Карло, красный - среднее значение всех вычисленных цен в
течение испытания.
22
Рис. 10: По горизонтали N, по вертикали ошибка вычислений по сравнению
со средним значением.
Как и в предыдущих случаях, графики показывают, что даже на небольшом количестве моделируемых испытаний ошибка метода уменьшается с
ростом количества случайных величин.
23
2
Параллелизм для ускорения вычислений
Так как при оценке опционов чрезвычайно важна скорость вычислений
(для быстрого реагирования на финансовом рынке), и в меньшей степени важна точность вычислений, метод Монте-Карло идеально подходит
для данной задачи. В методе Монте-Карло моделируется множество случайных величин, распределенных определенным образом, после чего находится среднее значение всех смоделированных величин, т.е. математическое ожидание. Таким образом, задачу моделирования случайных величин
можно разбить на части, которые могут выполняться параллельно, а затем сложить полученные результаты и найти среднее. Метод Монте-Карло
является задачей, относящейся к параллелизму по данным (SIMD - Single
Instruction Multiple Data), когда параллельные ядра графического ускорителя выполняют одну и ту же операцию над разными элементами данных.
В нашем случае, входные данные одинаковые, а разные только сгенерированные случайные величины.
Следовательно, целесообразно использовать технологии параллельных
вычислений для оценки стоимости опционов ради достижения ускорения
получения результата. Рассмотрим две технологии: автоматическое распараллеливание с помощью MATLAB Parallel Toolbox и технологию гибридных вычислений OpenCL.
2.1
Вычислительные мощности
Тестирование работы параллельных алгоритмов было проведено с использованием вычислительных ресурсов Ресурсного Центра "Вычислительный центр СПбГУ".[4] Вычисления в MATLAB проводились на кластере
Huawei, который обладает следующими характеристиками[16]:
∙ пиковая производительность 133 терафлопс;
∙ 20 ускорителей nVidia Tesla K40 (основанный на архитектуре NVIDIA
Kepler, 12 ГБ скоростной памяти DDR5, 2880 параллельных ядер)[17];
24
∙ коммутатор сети передачи сообщений с технологией InfiniBand, обеспечивающей скорость передачи данных в 56Гб/с.
Вычисления на OpenCL проводились на другом кластере, ввиду технических неполадок на первом. Это гибридный кластер, обладающий следующими характеристиками:
∙ 112 графических ускорителей NVIDIA Tesla M2050(448 параллельных
ядер, 3 ГБ DDR5);
∙ 2 порта Infiniband, обеспечивающих скорость передачи до 40Гб/с;
∙ всего 288 вычислительных ядер.
2.2
2.2.1
Параллельные вычисления в MATLAB
О MATLAB и Parallel Computing Toolbox
MATLAB - это высокоуровневый язык технических расчетов, а так же
интерактивная среда разработки алгоритмов и совокупность множества
пакетов прикладных программ для решения математических задач, разрабатываемая компанией MathWorks[18]. Особенно удобен MATLAB для
работы с матрицами и векторами, так как предоставляет основанные на
матрицах структуры данных. В данной работе были использованы следующие версии MATLAB: R2013b для последовательных алгоритмов и R2015b
для параллельных.
Возможность использовать для вычисления одной задачи несколько
ядер появилась в MATLAB благодаря пакету Parallel Computing Toolbox.
Он позволяет превратить алгоритм в параллельный гораздо проще, чем это
можно сделать с помощью других технологий. В версии MATLAB, которая
установлена на кластере, стоит ограничение на количество используемых
устройств, поэтому максимальное количество ядер, на которых возможно
провести расчеты, равно 32.
25
2.2.2
Параллельный алгоритм в MATLAB
Для работы с несколькими ядрами одновременно в MATLAB Parallel
Computing Tool существуют функции управления параллельным пулом.
Сначала нужно создать профиль кластера (локального или распределенного), заполнив все необходимые поля, с указанием максимального количества возможных ядер (NumWorkers), папки назначения, в которой будут
создаваться задачи, дополнительные командные строки для отправки задачи в очередь и другие параметры. После этого требуется пройти валидацию, то есть проверку кластера на доступность и выполнение задач. После
того как валидация пройдена, и профиль кластера установлен по умолчанию, можно начинать с ним работу.
Для создания пула параллельных потоков, используется функция parpool
(workers), где workers - количество используемых устройств. Этой функцией можно инициализировать переменную, назовем её poolobj. Тогда эта
переменная будет ассоциироваться с пулом потоков и по окончании вычислительных работ можно будет прервать работу параллельного пула функцией delete(poolobj).
Для создания параллельного алгоритма на основе последовательного в
MATLAB, необходимо использовать оператор параллельного цикла parfor,
который имеет столько итераций, на сколько ядер будет распределена задача, но выполняются они не последовательно, как в обычном цикле for, а
одновременно. Таким образом, для написания параллельных версий алгоритмов каждого из методов, требуется использовать цикл parfor с количеством итераций, равным количеству ядер, и внутренний цикл на 𝑁/𝑤𝑜𝑟𝑘𝑒𝑟𝑠
итераций, где 𝑁 - это количество моделируемых цен, а 𝑤𝑜𝑟𝑘𝑒𝑟𝑠 - число
ядер. В итоге, каждое ядро выполняет задачу по вычислению и сложению
𝑁/𝑤𝑜𝑟𝑘𝑒𝑟𝑠 цен. После этого остается только сложить данные с разных
ядер и усреднить, а так же дисконтировать значение.
Тогда формальный универсальный параллельный алгоритм можно представить следующим образом:
1. Инициализация переменных. Получение начальных значений, обну26
ление тех переменных, значения для которых появляются позже.
2. Инициализация объекта параллельного пула poolobj с количеством
требуемых потоков.
3. Запуск счетчика времени tic.
4. Цикл parfor i=1:workers.
5. Внутренний цикл for k = 1:N/workers, в теле которого соответствующим методом высчитывается сумма N/workers цен.
6. Конец цикла for, за ним конец цикла parfor.
7. "Свертка" - цикл for, складывающий все i-ые суммы с разных потоков. Можно так же реализовать функцией sum(A), если A - вектор.
8. Усреднение полученной суммы путем деления на N и дисконтирование полученный цены.
9. Остановка счетчика времени toc.
10. Завершение работы параллельного пула delete(poolobj).
2.3
2.3.1
Технология параллельных вычислений OpenCL
Описание технологии OpenCL
OpenCL - технология для гибридных вычислений, позволяющая проводить вычисления на различных процессорах (как CPU, так и GPU) параллельно. При этом процессоры могут так же быть и разных фирм, что
является достоинством OpenCL, так как нет ограничения в выборе аппаратуры: на данный момент появляются фирмы, изготавливающие, например,
графические ускорители, выдерживающие конкуренцию с самыми известными продуктами. Если технология CUDA поддерживает только графические процессоры Nvidia, то OpenCL может работать с ускорителями разного производства, а кроме них и с CPU (Central Processing Unit), и с другими
видами процессоров. OpenCL - относительно новая технология, разрабатывающаяся Khronos Group[19], но уже имеет несколько релизов.
27
В OpenCL структура вычислительной системы выглядит следующим
образом: есть один управляющий центральный процессор, называющийся
хостом (host), и набор процессоров, которые могут быть как ускорителями,
так и центральными процессорами, называющиеся вычислителями устройствами (computing devices). Хост устанавливает опции вычисления, организует работу и память устройств, отсылает программы для выполнения
на устройства, дожидается окончания работы, выполняет пост-обработку
данных и освобождает память. С другой стороны устройства получают
данные от хоста, выполняет программу, которая так же приходит от хоста,
и отсылает результат обратно.
В состав OpenCL входят собственный язык для записи кода, компилируемого и исполняемого на устройствах (так называемого kernel-кода или
device-кода), а так же API для записи кода для управляющего центрального процессора (host-кода). На данный момент API OpenCL поддерживает
разработку на языках C и C++, но так же существуют и сторонние API,
поддерживающие такие языки, как Java, Python, .NET и другие.
Что касается обмена данными, то в OpenCL есть иерархия памяти, которая состоит из следующих элементов[14]:
∙ Память хоста - память, доступная лишь управляющему процессору;
∙ Глобальная память - память, доступная как хосту, так и любому из
вычислительных устройств (определяется идентификатором __global);
∙ Постоянная память (константная) - часть глобальной памяти, менять
которую может только хост, устройства же только читать;
∙ Локальная память - область памяти, доступная только одному конкретному устройству(определяется идентификатором __local);
∙ Частная память (приватная) - память, доступная только одному вычислительному элементу устройства (определяется идентификатором
__privat).
Вычислительные элементы в OpenCL объединяются в группы, которые
называются Work group на уровне программы или Compute unit на уровне
28
аппаратуры. Каждый юнит и соответствующая группа имеет свою локальную память, которая гораздо быстрее общей, глобальной памяти, и поэтому
позволяет ускорить взаимодействие между потоками благодаря использованию элементами одной группы локальной, а не глобальной памяти для
передачи данных.
Уникально идентифицировать Work item можно либо по его глобальному порядковому номеру, либо по сочетанию координат группы в матрице
индексов и элемента в этой группе. Идентифицировать по одному элементу
в каждой группе можно используя локальный порядковый номер, который
начинается с 0 и до GroupSize-1.
2.3.2
Генератор случайных чисел
Как известно, для реализации метода Монте-Карло необходимо генерировать множество случайных величин. Однако в языке OpenCL C, который
используется для написания kernel’ов, не предусмотрен генератор случайных чисел. Хоть и основан на C, он не включает в себя такие библиотеки
как stdio и stdlib. В связи с этим в коде для вычислительных устройств
нельзя использовать даже самую простую функцию rand(), входящую в
состав языка C.
Кроме того, нет так же и других подключаемых библиотек для OpenCL,
которые бы предоставляли генератор случайных чисел. Таким образом,
приходится искать искусственные методы для выхода из данной ситуации.
Первый вариант - это генерировать весь массив случайных величин на хосте и пересылать исполнителям, а второй - использовать такие генераторы,
как Xorshift, Java random или Mersenne Twister, реализуя их самим.
Первый подход нежелателен для распараллеливания, когда требуется
как можно больше вычислений проводить на GPU независимо друг от друга. При реализации этого подхода, придется последовательно генерировать
весь набор случайных величин, которые понадобятся для расчета и записывать его в глобальную память. Кроме того, процессам придется каждый
29
раз обращаться в медленную глобальную память для считывания необходимого значения.
Для исследования влияния первого подхода на время работы программы, было реализовано динамическое выделение памяти для массива, так
как статического не хватало для записи уже 4 миллионов случайных величин. После этого, заводился цикл for для заполнения массива, в котором
вызывалась функция, генерирующая с помощью rand() и преобразования
Бокса-Мюллера стандартную нормально распределенную случайную величину, которая необходима для расчетов.
В проведенных экспериментах таким образом генерировалось до 40 миллионов таких случайных величин. Расчет времени выполнения программ
показал, что в итоге такой параллельный алгоритм дает ускорение только
в полтора раза, и примерно 90% времени занимает генерация массива. Например, при расчете азиатского опциона (СДУ) параллельным алгоритмом
затратилось время 118 секунд, а без учета генерации массива, 3 секунды.
Естественно, это доказывает, что использование такого метода возможно,
если количество требующихся случайных величин в пределах тысячи или
нескольких тысяч максимум. Но для метода Монте-Карло такой подход
бесполезен.
Поэтому, для решения этой задачи был реализован генератор с возможностью рассчитывать случайные величины прямо на GPU. Генераторы
случайных величин используют для своей активации некоторое значение
- "зерно"(seed), из которого затем различными операциями вычисляются
случайные величины. Прямо на GPU генерировать такое значение не представляется возможным: генерировать случайные числа там нельзя, если
взять за зерно время, может получиться, что у многих процессов одинаковое зерно, а если выразить значение с помощью номера процесса, каждый
раз будет генерироваться одно и то же. Поэтому необходимо генерировать
эти зерна на хосте.
Далее мы сталкиваемся с той проблемой, что если каждый процесс бу30
дет вычислять одну цену, то для каждого процесса придется отправлять
зерно и, соответственно, генерировать его последовательно на CPU. Тогда
размеры такого массива с инициализирующими значениями достигает масштабов массива, рассмотренного в первом подходе, и написание генератора
теряет смысл. Поэтому, было решено, что каждый процесс будет выполнять
𝑡𝑖𝑚𝑒𝑠 вычислений цен, тогда на хосте необходимо будет вычислять в 𝑡𝑖𝑚𝑒𝑠
раз меньше "зерен". Одного "зерна" будет достаточно для генерации необходимого количества случайных величин для одного процесса. Реализация
на практике показала, что такой метод является решающим проблему отсутствия генератора и больших временных затрат.
Конкретно, реализованный алгоритм основан на Java random[20]. С помощью него генерируются псевдослучайные равномерно распределенные
числа. Сам алгоритм берет seed типа unsigned long и производит с ним
некоторые преобразования, например, разного рода сдвиги, а затем вычленяет из него часть, которая станет очередным псевдослучайным числом.
Далее, как и в последовательной реализации, применяется алгоритм БоксаМюллера для преобразования в стандартное нормальное распределение.
Результаты замера ускорения параллельных программ с таким подходом к
генерации случайных величин рассмотрены в Главе 3.
2.3.3
Редукция выходного массива
При выполнении на множестве вычислительных элементов kernel-кода,
который использует метод Монте-Карло, на выходе получается большой
массив с данными, которые нужно обработать. Конечно, можно сделать
это на хосте, но быстрее будет провести редукцию при выполнении kernel’а,
когда потоки выполняются одновременно и имеют быстрый доступ к локальной памяти. Тогда число последовательных операций сократится.
Благодаря иерархии памяти и структуре организации системы в OpenCL,
реализовать данный подход представляется возможным. Так как все рабочие элементы группируются, и у каждой группы есть своя локальная память, то доступ к этой памяти элементы получают гораздо быстрее, чем
31
к глобальной. Тогда нужно использовать возможность обработать данные
в этой локальной памяти до того, как отправлять их в глобальную. Пример работы такой редукции, состоящей из локального и глобального этапа,
приведен на Рис.11[15].
Рис. 11: Схема работы метода MapReduce в OpenCL.
Следует учесть, что в языке OpenCL не поддерживаются динамические
структуры данных, поэтому использовать этот подход можно, только если
заранее известно количество данных, которые нужно будет обработать. В
нашем случае, каждый рабочий элемент на выходе будет иметь всего одно
число - сумму цен, соответственно количество информации, представляемой к редукции, будет равно количеству элементов в рабочей группе, а
затем количеству групп.
Таким образом, последовательность действий такова:
1. Каждый рабочий элемент выполняет расчет, а затем часть локальной
редукции - т.е. добавляет свой результат к общему результату группы.
2. Затем вызывается локальный барьер, который предотвращает продолжение выполнения программы до момента, когда все элементы в
группе не выполнят свою часть редукции.
32
3. В каждой группе рабочий элемент с локальным номером 0, добавляет
результат его группы к глобальной результирующей переменной.
4. Вызывается глобальный барьер до момента, когда все группы не закончат свою редукцию.
5. Если нужно, рабочий элемент с глобальным номером 0 генерирует
финальный результат из полученной суммы необходимым преобразованием.
Необходимо обратить внимание, что при такой редукции кроме использования барьеров требуется так же использовать атомарные операции, которые выполняются в качестве единичных транзакций. При вызове этой
операции одним потоком, блокируется область памяти, которая будет использоваться, и другие потоки не могут в это время что-либо сделать с этой
областью. Таким образом, не нарушится работа программы и не возникнет
проблем с одновременным использованием одной и той же информации.
В OpenCL такие операции реализованы с помощью функций, начинающихся со слова "atomic". Есть атомарные операции сложения, вычитания,
замены, а так же атомарный инкремент и декремент. В нашем случае было
бы удобно использовать атомарное сложение (atomic_add()), но к сожалению, эта операция реализована только для переменных целого типа, а наш
выходной результат будет типом с плавающей точкой. Поэтому приходится использовать операцию замены (atomic_xchg()), в качестве аргументов
которой указывается адрес на область памяти результирующей переменной
и значение суммы текущего значения этой переменной и добавляемого значения.
2.3.4
Структура программы на OpenCL
Программа на OpenCL должна обязательно включать в себя некоторую последовательность действий для подключения к системе и установки kernel’ов на выполнение. В первую очередь, необходимо определить
платформу, которая будет использоваться. Так как есть много реализаций OpenCL, некоторые из них работают с одними видами процессоров,
33
другие с другими, и на одной системе может быть несколько различных
платформ, требуется выбрать одну, которая работает с нужными устройствами. Это происходит посредством функции clGetPlatformIDs().
Следующий шаг - найти доступные устройства для запуска программы, удовлетворяющие заданным требованиям. Например, если требуются
устройства только GPU, то в аргументы функции clGetDeviceIDs() записать параметр cl_device_type_gpu. В случае cl_device_type_all будут
рассматриваться все возможные варианты.
Затем необходимо создать контекст и очередь команд функциями
clCreateContext() и clCreateCommandQueue() соответственно. Затем с
помощью функций clCreateProgramWithSource(), clBuildProgram() и
clCreateKernel() исходный kernel-код считывается в переменную типа
cl_program, затем строится исполняемый файл программы и создаётся вычислительный kernel.
После этого необходимо выделить память устройства для входных и выходных данных, расположенную в глобальной области. Затем заполнить её
входными данными, а так же установить значения дл аргументов функции kernel. Эти операции выполняются функциями clCreateBuffer() и
clSetKernelArg().
Когда все эти приготовления закончены, kernel ставится в очередь на
вычислительных элементах, а хост ждет окончания расчетов. Затем остается только считать выходные данные, провести с ними некоторые операции,
если требуется, и вывести.
Последняя стадия программы на OpenCL - освобождение памяти и завершение работы. Этот этап производится с помощью функций
clReleaseKernel(), clReleaseMemObject(), clReleaseProgram()
и других[21].
34
2.3.5
Параллельный алгоритм расчета опционов на OpenCL
Для вычисления стоимости опционов необходимо передать в kernel в качестве константных аргументов все используемые для расчета параметры:
стоимость базового актива в момент времени 𝑡0 , strike-цена, волатильность
и безрисковая процентная ставка, а так же время действия контракта. Кроме того, аргументами в kernel-функции так же должны быть указатели
на области памяти, подготовленные для входных и выходных данных. Создать переменную в локальной памяти для записи промежуточного результата редукции можно непосредственно в самой функции __kernel следующим образом: __local float local_result. Так как мы генерируем массив
случайных величин нормального распределения на хосте и передаем, то в
область глобальной памяти, выделенной для входных данных, необходимо записать этот массив посредством функции clEnqueueWriteBuffer().
Кроме значений для расчета цены базового актива и массива случайных
величин, необходимо так же передать количество цен, которое будет рассчитано каждым потоком.
Размер рабочей группы выбирается OpenCL автоматически под используемое вычислительное устройство. В нашем случае это 1024. Тогда размер
массива с входными данными будет 1024× кол-во групп × кол-во цен, рассчитанных каждым элементом, и таким же будет суммарное количество
рассчитанных цен для европейского опциона. Для азиатского опциона количество элементов массива будет в 𝑚 раз больше для того же количества
рассчитываемых цен, так как для расчета одной цены требуется провести
𝑚 промежуточных расчетов. То же самое касается метода, основанного на
континуальном интеграле.
Первым делом в kernel функции необходимо инициализировать переменную для результирующего значения. Это может сделать один из потоков, например, с глобальным номером 0. Глобальный номер можно определить с помощью функции get_global_id(). Затем определить локальную
переменную для результата редукции группы и так же инициализировать
её. Затем дождаться, когда все потоки группы проведут инициализацию и
начать вычисления. Для этого заводится цикл for от 𝑖 = 𝑙 до 𝑖 = 𝑙 + 𝑡𝑖𝑚𝑒𝑠 с
35
шагом 1 для европейского опциона методом стохастического дифференциального уравнения и азиатского методом уравнения в частных производных, где 𝑙 = 𝑔𝑒𝑡_𝑔𝑙𝑜𝑏𝑎𝑙_𝑖𝑑(0) * 𝑡𝑖𝑚𝑒𝑠, а 𝑡𝑖𝑚𝑒𝑠 - количество цен для расчета
каждым потоком. Для остальных рассматриваемых методов цикл будет от
𝑖 = 𝑙 до 𝑖 = 𝑙 + 𝑡𝑖𝑚𝑒𝑠 * 𝑚, но с шагом в m, а внутри будет еще один цикл
for от 𝑘 = 0 до 𝑘 = 𝑚 с шагом в 1. С такой организацией циклов в первом
случае для расчета цены берется 𝑖-ый элемент массива со входными данными, а для второго случая с двойным циклом (𝑖 + 𝑘)-ый. Формулы для
расчета цен подробно описаны в Главе 1.
В цикле мы рассчитываем сумму 𝑡𝑖𝑚𝑒𝑠 моделируемых цен, затем атомарной операцией добавляем среднее арифметическое этой суммы к локальному результату и дожидаемся завершения этой операции другими
членами локальной группы. Затем поток с локальным номером 0 (который
определяется функцией get_local_id()) добавляет локальный результат
редукции к выходной переменной, и после этого поток с глобальным номером 0 производит деление значения результирующей суммы на количество
групп. После этого хосту остается только принять данные и дисконтировать цену.
36
3
Анализ полученного ускорения вычислений
3.1
Ускорение в MATLAB
Для расчета ускорения для каждого метода был проведен ряд испытаний, в каждом из которых моделировалось 1.200.000 случайных величин
для наиболее точного значения цены опциона и использовалось различное
количество ядер - от 1 до 32. Время работы замерялось с помощью функций tic (в момент начала отсчета) и toc(в момент окончания).
Благодаря MATLAB, было получено значительное ускорение для каждого из методов, не смотря на то, что количество ядер было сильно ограничено. В результате расчетов скорость выполнения программы увеличилась
почти в 100 раз для расчета азиатского опциона методом континуального
интеграла на 32 ядрах по сравнению с последовательным алгоритмом(см.
Рис.12).
Рис. 12: График полученного ускорения для расчет азиатского опциона
методом континуального интеграла в MATLAB.
37
Заметим, что ускорение для алгоритма расчета европейского опциона
этим же методом не было рассчитано, так как в них требуется провести
одинаковое количество расчетов, а следовательно и ускорение будет примерно одинаковое. Таким образом, для европейского опциона результаты
аналогичны приведенному графику для азиатского(Рис.12).
Рис. 13: График полученного ускорения для расчет азиатского опциона
методом СДУ в MATLAB.
Для расчета цены азиатского опциона с помощью метода стохастического дифференциального уравнения было получено ускорение в 65 раз
на 32 потоках(Рис.13), а для европейского опциона ускорение достигло
92(Рис.14).
Для метода основанного на разностной схеме и блужданию по сетке,
вычисление стоимости азиатского опциона с помощью параллельного алгоритма на 32 ядрах заняло в 86 раз меньше времени, чем аналогичные
последовательные вычисления(Рис.15).
38
Рис. 14: График полученного ускорения для расчета европейского опциона
методом СДУ в MATLAB.
Рис. 15: График полученного ускорения для расчета азиатского опциона
методом уравнения в частных производных в MATLAB.
39
3.2
Ускорение при помощи OpenCL
Последовательные алгоритмы для сравнения с параллельным на OpenCL,
были написаны на языке C. Для генерации случайных величин в последовательном алгоритме использовалась функция rand(), предоставляемая библиотекой stdlib.h и преобразование Бокса-Мюллера, генерирующее из двух
равномерно распределенных случайных величин одну нормально распределенную. В остальном реализация алгоритмов, описанных в разделе 1.2,
на языке С не вызывает сложностей. Процесс реализации параллельного
алгоритма на OpenCL был достаточно подробно описан во второй главе.
Для замера времени исполнения программ, как последовательных, так
и параллельных, была использована функция gettimeofday() стандартной библиотеки time.h и некоторые преобразования с типом для вывода в
удобном виде.
Европейский опцион(СДУ)
Рис. 16: График полученного ускорения для расчета европейского опциона
в OpenCL(СДУ).
40
Рассмотрим ускорение, полученное при разном количестве усредняемых
генерируемых цен, для алгоритма оценки европейского опциона первым из
рассматриваемых подходов (стохастическое дифференциальное уравнение)
на графике, представленном на Рис.15. На нем изображено ускорение параллельного алгоритма относительно последовательного.
Как видно на графике, ускорение растет с ростом числа случайных
величин и достигает значения в 80 при 40 миллионах. Когда проводились
опыты с генерацией случайных величин на CPU, было выяснено, что с учетом времени генерации ускорение составляло всего 1,5, но даже без учета
времени генерации, ускорение достигло только 60. Этот пример показывает,
как важно использовать локальную и частную память, так как обращаться за данными в глобальную область памяти занимает много времени и
отрицательно сказывается на ускорении. При сравнении ускорения с ускорением MATLAB стоит помнить, что по техническим причинам, задачи на
OpenCL считались на более слабых GPU, и нужно делать на это скидку.
Рассмотрим теперь ускорение, полученное при использовании 1, 2 и 3
плат GPU, при расчете той же задачи(Рис.16).
Рис. 17: График полученного ускорения для расчета европейского опциона
в OpenCL(СДУ) на разном количестве GPU.
41
Так, если на одном GPU ускорение на больших количествах расчетов
равно примерно 80, то на двух уже примерно 110, а на трёх достигает
150(Рис.17). Такая хорошая масштабируемость является результатом того,
что вычисления на разных поток почти не зависят друг от друга(процессы
синхронизируются только во время редукции). При этом, результаты для
двух и трех GPU могли бы быть еще лучше, если бы нам не приходилось
для них последовательно рассчитывать, соответственно, в два и в три раза
больше зерен, чем для 1 GPU.
Азиатский опцион(СДУ)
Теперь рассмотрим аналогичные графики для азиатского опциона (СДУ).
Для его расчета, в отличие от европейского, необходимо дополнительно
вычислять 𝑚 промежуточных сумм, а значит, каждый поток делает в 𝑚
раз больше вычислений. Поэтому время вычисления сразу увеличивается
в несколько раз, а ускорение не достигает таких больших значений на соответственных количествах случайных величин.
На одном GPU ускорение почти всегда держится на уровне 56 раз(Рис.18).
Для двух и трёх GPU ускорение не сильно возросло по сравнению с одним(Рис.19). Это связано с тем, что для большего количество GPU нужно
генерировать больше значений зерен на CPU.
Европейский опцион(КИ)
Рассмотрим теперь ускорение, полученное при проверке работы параллельной реализации алгоритма расчета европейского опциона через континуальный интеграл. Как и в разделе про ускорение в MATLAB, мы будем
рассматривать ускорение оценки стоимости только одного вида опционов
из рассматриваемых нами, так как алгоритмы для них отличаются всего
одной операцией и требуют одинаковое количество вычислений, а это означает, что время работы будет схоже или даже одинаково.
42
Рис. 18: График полученного ускорения для расчета азиатского опциона в
OpenCL(СДУ).
Рис. 19: График полученного ускорения для расчета азиатского опциона в
OpenCL(СДУ) на разном количестве GPU.
В результате замера времени, программа, реализованная на OpenCL и
работающая на 1 GPU, показала ускорение в 60-62 раза(Рис.20). Экспери43
менты так же показали, что и масштабируемость у этого алгоритма тоже
неплохая - на трех GPU удалось достигнуть ускорение в 90 раз(Рис.21).
Рис. 20: График полученного ускорения для расчета европейского опциона
в OpenCL(КИ).
Эти результаты сходны с результатами, полученными для азиатского
опциона, считаемого через стохастическое дифференциальное уравнение.
Этот факт объясняется тем, что алгоритмы имеют похожую структуру,
разница заключается лишь в способе расчета цен, а количество вычислений одинаковое.
44
Рис. 21: График полученного ускорения для расчета европейского опциона
в OpenCL(КИ) на разном количестве GPU.
Азиатский опцион(УвЧП)
Рис. 22: График полученного ускорения для расчета азиатского опциона в
OpenCL(УвЧП).
45
Данный метод показал себя хорошо в плане ускорения - почти такое
же ускорение, как для европейского опциона через СДУ. На одном GPU
ускорение составило 70 (Рис.22), а на трех GPU достигло почти 130(Рис.23).
Рис. 23: График полученного ускорения для расчета азиатского опциона в
OpenCL(УвЧП) на разном количестве GPU.
3.3
Сравнение времени работы реализаций в MATLAB,
на языке C и OpenCL
В предыдущих разделах мы рассмотрели ускорение алгоритмов в
MATLAB и OpenCL и заметили, что ускорение, полученное в OpenCL,
немного ниже, чем в MATLAB. Это объясняется отсутствием в языке
OpenCL генератора псевдослучайных величин и , следовательно, необходимостью генерировать зерна для написанного генератора на CPU последовательно, а так же тем, что по техническим причинам вычисления проводились на разных кластерах. MATLAB работал на более мощном кластере.
Однако, стоит заметить, что в рассматриваемой задаче важно не только получить ускорение, но и следить за временем выполнения программы.
46
Так, что последовательная, что параллельная, реализация алгоритмов в
MATLAB работает гораздо медленнее реализаций тех же алгоритмов на C
и OpenCL. Это связано с более высоким уровнем организации MATLAB, и,
следовательно большим количество надстроек, которые работают дольше.
Языки C и OpenCL являются более низкоуровневыми.
Например, при 1200000 случайных величин последовательная программа для европейского опциона в MATLAB работает 130 секунд, в то время
как реализация на языке C на таком же количестве случайных величин
работает около 0,4с. В таких же условиях последовательная программа
для азиатского опциона на MATLAB работает 176 секунд, а на C около
5 секунд. При этом, самое быстрое параллельное выполнение первой программы на MATLAB 1.4с, а второй 2.7с. При параллельной реализации на
OpenCL на 1 GPU европейский опцион считается за 0,005 секунд, а азиатский за 0,085c.
47
Выводы
В результате проведения исследования методов и анализа полученных
данных были сделаны следующие выводы:
1. Время работы реализаций в MATLAB больше, чем на языке Си и
OpenCL. Это объясняется более высоким уровнем организации, а так
же тем, что язык MATLAB является интерпретируемым.
2. Однако ускорение работы параллельного алгоритма относительно последовательного в MATLAB достаточно высоко. При 32-х потоках было достигнуто ускорение почти в 100 раз для двух из трех методов.
Можно сделать вывод, что параллельные алгоритмы в MATLAB работают хорошо, и удобно использовать MATLAB Parallel Computing
Toolbox, если алгоритм требует работы с матрицами или автор алгоритма не имеет достаточных знаний в языках программирования.
3. Последовательные реализации алгоритмов на языке C гораздо быстрее аналогичных на MATLAB, поэтому в условиях рассматриваемой
задачи целесообразно при возможности использовать реализации на
С или других подобных языках, так как главной целью написания
программ, оценивающих стоимость опционов, является как можно
более быстрое получение результата. Кроме того, в рассматриваемых
алгоритмах не требуется операций с матрицами или других сложных
математических вычислений, поэтому нет причин выбирать реализацию на MATLAB.
4. Очень большую роль в скорости выполнения программы играет генерация случайных величин, так как их требуется огромное количество. Поэтому, программы на OpenCL выполняются немного медленнее, чем могли бы, так как приходится генерировать массив "зерен"размером в количество используемых процессов - для каждого
процесса одно "зерно". Этот факт замедляет работу и уменьшает
ускорение относительно последовательного алгоритма, но ускорение
всё равно есть и заметное.
5. OpenCL достигает ускорения в 80 раз на одной плате GPU и до 150
48
раз на трёх платах GPU при оценке европейского опциона (СДУ).
Для азиатского опциона(СДУ) ускорение вышло 56 раз на одной плате и 72 раза на трех платах GPU. Для алгоритма расчета через континуальный интеграл для обоих видов опционов ускорение на OpenCL
достигло 62 раз на одной плате и 90 раз на трех платах GPU. И для
расчета азиатского опциона через уравнение в частных производных,
ускорение достигло 70 на одной плате GPU и 130 на трёх платах GPU.
49
Заключение
В процессе рассмотрения задачи, были разобраны и реализованы три
подхода к нахождению цены европейского и азиатского Call-опциона (таким же образом можно искать и Put-опцион, формула его стоимости немногим отличается от Call, главное - правильно спрогнозировать цену на базовый актив в нужный нам момент времени). Все эти методы сходятся и
дают нам примерно один и тот же результат достаточной точности.
Для получения ускорения работы алгоритмов были написаны программы на MATLAB, C и OpenCL, которые запускались на высокопроизводительных кластерах. В результате измерения времени данных алгоритмов
было получено различное ускорение.
Реализации на MATLAB дали хорошее ускорение, но время работы программ было неудовлетворительным. Реализации на OpenCL показали чуть
меньшее ускорение ввиду особенностей работы написанного для алгоритмов генератора, но время работы программ оказалось гораздо меньше.
В итоге были выполнены все поставленные задачи и сделан вывод, что
ускорение расчета опционов на OpenCL имеет потенциал.
50
Список литературы
[1] Ермаков С. М.
Метод Монте-Карло в вычислительной математике:
Вводный курс. СПб.: Невский Диалект; М.: БИНОМ. Лаборатория
знаний, 2009. 192 c.
[2] Войтишек А. В. Михайлов Г.А.
Численное статистическое модели-
рование. Методы Монте-Карло. Академия М., 2006.
[3] Дегтярев А.Б. Андрианов С.Н.
Параллельные и распределенные вы-
числения. Часть 1. Спб.:"СОЛО" 2007. 60 с.
[4] Сайт Ресурсного Центра "Вычислительный центр СПбГУ". http://
cc.spbu.ru/.
[5] Fischer Black and Myron Scholes. The pricing of options and corporate
liabilities.
The journal of political economy, pages 637–654, 1973.
[6] Hongbin Zhang.
Pricing Asian Options using Monte Carlo Methods.
Department of Mathematics Uppsala University, 2009. 36 c.
[7] Михаил
Глухов.
Оценка
опционов
методом
Монте-Карло.
Futures&Options, апрель 2009. 38-43 с.
[8] Михаил Глухов.
Карло.
Оценка экзотических опционов методом Монте-
Futures&Options, май 2009. 40-49 с.
[9] Guido Montagna, Oreste Nicrosini, and Nicola Moreni. A path integral way
to option pricing.
Physica A: Statistical Mechanics and its Applications,
310(3):450–466, 2002.
[10] Vadim Linetsky. The path integral approach to financial modeling and
options pricing.
Computational Economics, 11(1-2):129–163, 1997.
[11] Daniel Zwillinger.
Handbook of differential equations, volume 1. Gulf
Professional Publishing, 1998.
Опционы, фьючерсы и другие производные финансовые
инструменты, 6-е издание. Издательский дом Вильямс, 2008.
[12] Джон К Халл.
51
[13] Центр компетенций mathworks matlab.exponenta.
http:\matlab.
exponenta.ru.
[14] Aaftab Munshi, Benedict Gaster, Timothy G Mattson, and Dan Ginsburg.
OpenCL programming guide. Pearson Education, 2011.
[15] Matthew Scarpino. Opencl in action: How to accelerate graphics and
computation. ny, 2012.
[16] Сайт Научного Парка СПбГУ. http://researchpark.spbu.ru.
[17] Сайт компании nvidia. http://www.nvidia.ru.
[18] Сайт the mathworks. www.mathworks.com.
[19] Сайт компании khronos group. https://www.khronos.org/opencl/.
[20] Электронная документация oracle по java. http://docs.oracle.com/
javase/7/docs/api/index.html.
[21] Электронная документация khronos group по opencl. https://www.
khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/.
52
Отзывы:
Авторизуйтесь, чтобы оставить отзыв