Санкт-Петербургский государственный университет
Кафедра компьютерного моделирования и многопроцессорных систем
Попова Евгения Николаевна
Выпускная квалификационная работа бакалавра
Численный анализ структур в финансовой математике
Направление 010300
Фундаментальная информатика и информационные технологии
Научный руководитель,
доктор физ.-мат. наук,
профессор
Богданов А. В.
Санкт-Петербург
2016
Содержание
Введение ......................................................................................................................................................3
Постановка задачи ......................................................................................................................................5
Обзор литературы .......................................................................................................................................9
Глава 1. Модель Блэка-Шоулза. ...............................................................................................................10
Глава 2. Методы Монте-Карло ................................................................................................................12
2. 1. Стохастическое дифференциальное уравнение .........................................................................12
2. 1. 1. Стохастическое дифференциальное уравнение для европейского опциона ..................12
2. 1. 2. Стохастическое дифференциальное уравнение для азиатского опциона ......................13
2. 2. Континуальный интеграл ..............................................................................................................14
2. 2. 1. Континуальный интеграл для европейского опциона .......................................................14
2. 2. 2. Континуальный интеграл для азиатского опциона ............................................................16
Глава 3. Обзор вычислительных систем ..................................................................................................18
3. 1. Обзор Nvidia CUDA .........................................................................................................................19
Глава 4. Реализация с помощью CUDA ....................................................................................................22
4. 1. Алгоритм параллельной редукции ..............................................................................................22
4. 2. Алгоритм префиксных сумм .........................................................................................................22
4. 3. Генерация случайных чисел. ........................................................................................................24
Глава 5. Численные эксперименты. .........................................................................................................26
5.1. Вычисление цены европейского опциона с помощью стохастического дифференциального
уравнения ...............................................................................................................................................27
5. 2. Вычисление цены азиатского опциона с помощью стохастического дифференциального
уравнения ...............................................................................................................................................30
5. 3. Вычисление цены европейского опциона с помощью континуального интеграла ................32
5. 4. Вычисление цены азиатского опциона с помощью континуального интеграла .....................33
5. 5. Масштабируемость .......................................................................................................................34
Выводы .......................................................................................................................................................37
Заключение ................................................................................................................................................38
Список литературы ....................................................................................................................................39
2
Введение
Ценообразование опционов является одной из важных задач в
финансовой математике. В настоящее время объем вычислений в различных
задачах становится все больше, и архитектура вычислительных систем
постоянно
меняется,
поэтому
эффективная
реализация
моделей
ценообразования опционов на современных устройствах становится все
более важной проблемой. Более того, огромное количество задач со
сложными математическими моделями не могут быть решены аналитически.
В таких случаях используются численные методы, которые требуют высокой
вычислительной мощности. Чтобы получить достаточно точные результаты,
временные
затраты
могут
быть
большими,
поэтому
ускорение
на
графических процессорах является эффективным способом решения задач с
помощью численных методов.
Технология CUDA, разработанная NVIDIA, является архитектурой
параллельных
вычислений
общего
назначения,
которая
использует
графические процессоры для решения сложных задач. Архитектура CUDA
позволяет разработчикам использовать языки C и C++, которые являются
наиболее широко используемыми языками программирования высокого
уровня.
В финансовой сфере время играет огромную роль, так как любая
задержка в обработке информации может привести к экономическим
потерям. Поэтому использование параллельных вычислений при построении
модели ценообразования опционов является оправданным.
Огромную роль в открытии методов вычисления цены европейских
опционов сыграла работа [1]. Формула Блэка-Шоулза является одним из
основных вычислительных инструментов, который измеряет и прибыль, и
убытки, и риск опционной сделки для инвесторов. Для уравнения Блэка3
Шоулза
можно
найти
точное
аналитическое
решение,
однако
для
обобщенной модели в некоторых случаях (напр. в случае азиатского
опциона) точное решение найти невозможно. Эта проблема приводит к
использованию численных методов, которые в свою очередь приводят к
приближенному значению цены опциона.
Часто используют методы Монте-Карло [2][3], которые основаны на
получении большого числа стохастических процессов. В настоящее время
финансовые расчеты часто основываются на этих методах из-за присущей им
высокой степени параллелизма.
Таким образом, исследование оценки опционов с помощью методов
Монте-Карло делают данную работу актуальной.
4
Постановка задачи
Для
дальнейшего
понимания
работы
стоит
дать
следующие
определения [4]:
Опр. 1. Опцион – это контракт, согласно которому покупатель опциона
получает право (не обязательство) купить или продать базовый актив по
заранее оговоренной цене, которая называется ценой «страйк» или ценой
исполнения, в определенный момент в будущем или на протяжении
определенного отрезка времени. В то же время, продавец опциона обязан
совершить сделку по продаже или купле базового актива в соответствии с
условиями проданного опциона.
Опр. 2. «Колл»-опцион – опцион, предоставляющий покупателю право
совершить покупку базового актива по цене «страйк».
Опр. 3. «Пут»-опцион – опцион, предоставляющий покупателю право
продать базовый актив по цене «страйк».
Условия
исполнения
опциона
(продажи
или
покупки
актива)
принимается продавцом опциона, а покупатель должен соблюдать эти
условия. Следовательно, корректное значение опциона имеет решающее
значение для совершения сделки для обеих сторон.
Введем следующие обозначения:
𝐶(𝑆, 𝑡) – цена опциона, зависящая от времени и цены базового актива,
𝑆(𝑆0 ) – цена базового актива (начальная цена базового актива),
𝜎 – волатильность базового актива, статистический финансовый показатель,
который характеризует изменчивость цены базового актива,
𝑟 – безрисковая процентная ставка,
5
𝐴(𝑆, 𝑡) – средняя цена опциона за определенный временной промежуток,
зависящая от времени и цены базового актива,
𝐾 – цена-«страйк» опциона,
𝑇 – время экспирации (момент или временной отрезок, в который может
быть выполнен опцион),
𝑡 – текущее время.
Существуют различные виды опционов [5]. Наиболее популярными
являются американский, европейский и некоторые виды экзотических
опционов (азиатский, бинарный, барьерный и др.). Рассмотрим некоторые из
них подробнее:
1. Европейский опцион может быть исполнен только в момент
экспирации. Тогда в случае колл-опциона, если цена базового актива 𝑆 к
моменту истечения контракта будет выше, чем цена исполнения 𝐾, тогда
покупатель, при совершении сделки, получит прибыль, равную 𝑆𝑇 − 𝐾. Если
же цена базового актива в момент экспирации будет ниже цены исполнения,
то покупатель не получит никакой прибыли, так как опцион будет не
выгодно исполнять.
Тогда в этом случае справедливая цена опциона, также называемая
премией,
которая
выплачивается
продавцу,
следующим образом:
В случае европейского колл-опциона:
𝐶 ( 𝑇) = max(𝑆𝑇 − 𝐾, 0).
В случае европейского пут-опциона:
𝐶( 𝑇) = max(𝐾 − 𝑆𝑇 , 0).
6
может
быть
вычислена
2. Азиатский опцион – это опцион, цена исполнения которого
основывается на средней стоимости базового актива за определенный
временной промежуток. Преимуществом этого опциона является то, что
продавцу
намного
сложнее
управлять
условиями
контракта,
чтобы
уменьшить прибыль покупателя, так как она будет зависеть только от цены
базового актива в моменты мониторинга и не будет фиксированной и заранее
известной в момент покупки опциона, как в случае с европейским типом.
Тогда справедливая цена азиатского опциона может быть вычислена
следующим образом:
1) Если мониторинг цены базового актива является непрерывным:
Для колл-опциона
𝑇
1
𝐶( 𝑇) = max( ∫ 𝜔(𝑡)𝑑𝑡 − 𝐾 ,0).
𝑇
0
Для пут-опциона
𝑇
1
𝐶( 𝑇) = max(𝐾 − ∫ 𝜔(𝑡)𝑑𝑡 ,0).
𝑇
0
2) Если мониторинг цены базового актива является дискретным:
Для колл-опциона
𝑛
1
𝐶( 𝑇) = max( ∑ 𝜔(𝑡𝑖 ) − 𝐾 ,0).
𝑇
𝑖=0
Для пут-опциона
𝑛
1
𝐶( 𝑇) = max (𝐾 − ∑ 𝜔(𝑡𝑖 ) ,0) .
𝑇
𝑖=0
Существуют различные модели ценообразования опционов:
7
Модель Блэка-Шоулза
Биномиальная модель
Модель Хестона
Модель Монте-Карло
Модель Кокса-Рубинштейна
Другие
В этой работе будет использован метод Монте-Карло, а алгоритмы
нахождения цены будут реализованы с помощью:
стохастического дифференциального уравнения,
континуального интеграла.
Метод Монте-Карло
основан
на получении большого числа
стохастических процессов. Этому методу присуща высокая степень
параллелизма.
Поэтому
целью
ценообразования
этой
опционов,
работы
используя
является
построить
гетерогенные
модель
вычислительные
системы, и за счет этого снизить вычислительные издержки методов оценки
опционов.
Чтобы достичь положительного результата, необходимо выполнить
следующие задачи:
1. Исследовать алгоритмы методов оценки опционов;
2. Реализовать последовательный и параллельный алгоритмы, используя
CUDA;
3. Провести численные эксперименты и проанализировать их результаты;
8
Обзор литературы
В источниках [1], [4], [5] описаны основы финансовой математики, в
том числе связанные с опционами. Рассмотрены различные модели
ценообразования и фондовый рынок в целом.
В работах [2], [3] представлены описания методов Монте-Карло, в [6]
методы для стохастического дифференциального уравнения для азиатских и
европейских опционов.
В статьях [7], [8] исследованы многомерные интегралы, используемые
при ценообразовании опционов, которые находят применение в вычислении
цены азиатских опционов.
Работы
[11],
[12],
[14]
помогают
в
понимании
архитектуры
вычислительных систем и параллельных вычислений и в освоении
технологии CUDA.
Источники [13], [15], [16] являются официальными и предоставлены
компанией NVIDIA.
9
Глава 1. Модель Блэка-Шоулза.
Перед тем как решить любое уравнение, стоит попытаться решить его
аналитически. И только в том случае, если аналитического решения не
существует или его решение является трудоемкой задачей, используются
численные методы, которые сводятся к выполнению конечного числа
действий над числами. В таком случае полученный результат имеет
определенную погрешность, в то время как аналитическое решение является
точным.
Уравнение Блэка-Шоулза позволяет получить аналитическое решение.
Впервые модель Блэка-Шоулза для оценки опционов была опубликована в
статье [1] в 1973 году. Эта формула была открыта экономистами Фишером
Блэком, Мирном Шоулзом и Робертом Мертоном. Данную модель
используют и в настоящее время.
Ключевым элементом при оценке опциона с помощью этой модели
является ожидаемая волатильность базового актива. Так как возрастающая
или убывающая цена на него влияет прямо пропорционально на цену
опциона.
Однако для использования данной модели должны быть выполнены
следующие условия:
В течение всего срока действия опциона не выплачиваются
дивиденды по базовому активу.
Опцион являются европейским, то есть он может быть исполнен
только в момент экспирации.
Финансовые рынки являются эффективными.
Комиссий и транзакционных издержек нет.
Безрисковая процентная ставка и волатильность являются
постоянными и заранее известными.
10
Модель основана на логнормальном распределении базового
актива.
Тогда формула Блэка – Шоулза выглядит следующим образом:
𝐶(𝑆, 𝑡) = 𝑆𝑁(𝑑1 ) − 𝐾𝑒 −𝑟(𝑇−𝑡) 𝑁(𝑑2 )
𝑆
𝜎2
ln ( ) + (𝑟 + ) (𝑇 − 𝑡)
𝐾
2
𝑑1 =
𝜎 √𝑇 − 𝑡
(1)
𝑑2 = 𝑑1 − 𝜎√𝑇 − 𝑡
𝑆𝑁(𝑑1 ) характеризует ожидаемую прибыль от покупки базового актива.
𝐾𝑒 −𝑟(𝑇−𝑡) 𝑁(𝑑2 ) определяет размер цены исполнения в момент экспирации.
Необходимо
заметить,
что
допущения,
изложенные
выше,
соответствует идеальному рынку, а значит, их полное выполнение не
является реальным. Однако именно эта модель послужила началом для
дальнейшего исследования опционов на фондовом рынке.
11
Глава 2. Методы Монте-Карло
2. 1. Стохастическое дифференциальное уравнение
Опр.
4.
Стохастическое
дифференциальное
уравнение
(СДУ)
—
дифференциальное уравнение, в котором один член или более представляют
собой стохастический процесс (случайный процесс).
Соответствующее
модели
Блэка-Шоулза
(1)
стохастическое
дифференциальное уравнение выглядит следующим образом [5]:
𝑑𝑆(𝑡) = 𝑟𝑆(𝑡)𝑑𝑡 + 𝜎𝑑𝐵(𝑡),
где 𝐵(𝑡) – стандартный винеровский процесс, 𝐵(𝑡)~𝒩(0, 1).
Тогда цену базового актива в момент экспирации можно вычислить по
следующей формуле:
𝑆(𝑇) = 𝑆0 𝑒
𝜎2
)𝑇+𝜎 √𝑇𝜀
2
(𝑟−
(1)
где 𝜀 – случайная величина, распределенная по нормальному закону.
2. 1. 1. Стохастическое дифференциальное уравнение для
европейского опциона
Для европейских опционов метод Монте-Карло можно описать
следующим алгоритмом [6]:
1. Сгенерировать 𝑛 случайных чисел для 𝜀;
2. Для каждого значения 𝜀, вычислить цену базового актива 𝑆(𝑇),
по формуле (2)
3. Найти среднее значение цены базового актива 𝑆ср
4. Найти цену опциона по следующей формуле:
12
𝐶 = 𝑒 −𝑟𝑇 𝑆ср
2. 1. 2. Стохастическое дифференциальное уравнение для
азиатского опциона
Сложность применения метода Монте-Карло к ценообразованию
азиатских опционов заключается в том, что необходимо в каждый момент
наблюдения генерировать случайный временной ряд для цены базового
актива, в отличие от европейских опционов, где было достаточно
сгенерировать случайный ряд для цены базового актива только в момент
экспирации.
Для начала, разобьем время мониторинга на 𝑚 отрезков, то есть шаг
𝑇
∆𝑡 = :
𝑚
𝑡𝑖 =
1.
𝑖𝑇
,
𝑚
𝑖 = 0, 1, … , 𝑚
генерируем 𝑚 случайных значений 𝜀, и каждому числу
сопоставляем значение цены базового актива по следующей формуле:
𝑆(𝑡𝑖+1 ) = 𝑆(𝑡𝑖 )𝑒
2.
(𝑟−
𝜎2
) ∆𝑡+𝜎√ ∆𝑡𝜀
2
,
𝑖 = 0, 1, … , 𝑚
Находим среднее арифметическое значение полученного ряда:
𝑚
1
𝑆̅ = ∑ 𝑆(𝑡𝑖 )
𝑚
𝑖=1
3.
Повторяем 1 и 2 шаги 𝑛 раз. Получаем 𝑛 значений цен базового
актива.
4.
Вычисляем цену опциона:
𝑛
1
𝐶(𝑆0 , 𝑇) = 𝑒 −𝑟𝑇 ∑ 𝑆𝑗
𝑛
𝑗=1
13
2. 2. Континуальный интеграл
Опр. 4. Континуальный интеграл/функциональный интеграл – результат
функционального интегрирования, то есть вычисление интеграла некоторого
функционала Φ по пространству функций 𝑥(𝑡).
2. 2. 1. Континуальный интеграл для европейского опциона
Уравнение Блэка-Шоулза в частных производных имеет вид:
𝜕𝐶 1 2 2 𝜕 2 𝐶
𝜕𝐶
+ 𝜎 𝑆
+
𝑟𝑆
− 𝑟𝑉 = 0
𝜕𝑡 2
𝜕𝑆 2
𝜕𝑆
Сделаем следующую замену:
𝑥 = ln 𝑆
Тогда оно примет следующий вид [7]:
𝜕𝐶 1 2 𝜕 2 𝐶
𝜕𝐶
+ 𝜎
+
𝜇
− 𝑟𝐶 = 0
𝜕𝑡 2 𝜕𝑥 2
𝜕𝑥
𝜎2
𝜇=𝑟−
2
𝐶(𝑒 𝑋𝑇 , 𝑇) = 𝐹(𝑆𝑇 ) = max(𝑒 𝑋𝑇 − 𝐾, 0)
Для решения этой задачи можно использовать формулу Фейнмана-Каца:
𝐶(𝑆, 𝑡) = 𝑒 −𝑟𝑇 𝐸[𝐹(𝑆𝑇 )],
𝐸[𝐹(𝑆𝑇 )] – мат. ожидание, которое можно представить также в виде:
+∞
𝐸[𝐹(𝑆𝑇 )] = ∫−∞ max{ 𝑒 𝑥𝑇 − 𝐾, 0}𝑝(𝑥𝑇 , 𝑇|𝑥, 𝑡)𝑑𝑥𝑇 ,
где 𝑝(∙) – плотность вероятности перехода из (𝑥, 𝑡) в (𝑥𝑇 , 𝑇).
14
(3)
𝑝(𝑥𝑇 , 𝑇|𝑥, 𝑡) =
+∞
+∞
= ∫ ⋯ ∫ 𝑑𝑥1 ⋯ 𝑑𝑥𝑛
−∞
−∞
1
√(2𝜋𝜎 2 ∆𝑡)𝑛+1
𝑒
1
(− 2 ∑𝑛+1
(𝑥 −(𝑥𝑘−1 +𝜇∆𝑡))2
2𝜎 ∆𝑡 𝑘=1 𝑘
Выражение (3) можно представить в виде [8]:
+∞
+∞
𝐸[𝐹(𝑆𝑇 )] = ∫ ⋯ ⋯ ∫ 𝑑𝑥1 ⋯ 𝑑𝑥𝑛+1 max{𝑒 𝑥𝑛+1 − 𝐾, 0} ∙
−∞
−∞
1
𝑒
√(2𝜋𝜎 2 ∆𝑡)
(−
1
∑𝑛+1(𝑥 −(𝑥𝑘−1 +𝜇∆𝑡))2
2𝜎 2 ∆𝑡 𝑘=1 𝑘
𝑛+1
2
Тогда алгоритм вычисления цены европейского опциона методом
Монте-Карло с помощью континуального интеграла будет выглядеть
следующим образом:
1.
Вычислить 𝑥0 = ln 𝑆0 .
2.
Сгенерировать значения 𝑥1 , 𝑥2 , … , 𝑥𝑛+1 с помощью значения 𝑥0 и
генеральной совокупности с плотностью:
1
√2𝜋𝜎 2 ∆𝑡
𝑒
(−
1
(𝑥 −(𝑥𝑖−1 +𝜇∆𝑡))2
2𝜎 2 ∆𝑡 𝑖
3. Вычислить max{𝑒 𝑥𝑛+1 − 𝐾, 0}.
4. Повторить 1 – 3 шаги 𝑁 раз.
5. Полученные на третьем шаге значения обозначить как 𝑎𝑖 , 𝑖 =
1, … , 𝑁.
6. Найти значение цены опциона:
15
𝑁
𝑒 −𝑟𝑇
𝐶(𝑆0 , 𝑇) =
∑ 𝑎𝑖
𝑁
𝑖=1
7. Повторять все шаги до тех пор, пока значение, полученное на текущем
шаге, 𝐶(𝑆0 , 𝑇) не будет отличаться от значения, полученного на предыдущем
шаге, на достаточно малую величину.
2. 2. 2. Континуальный интеграл для азиатского опциона
Цена азиатского колл-опциона равна следующему выражению:
𝐹(𝑆, 𝐴) = max(𝐴(0, 𝑇) − 𝐾 ,0).
Согласно материалам из [] это значение можно также вычислить с
помощью следующей формулы:
𝐶(𝑆𝑇 , 𝑇) = 𝑒 −𝑟𝑇 𝐸[𝐹(𝑒 𝑥𝑇 , 𝐴(0, 𝑇))],
+∞
+∞
𝐸[𝐹(𝑒 𝑥𝑇 , 𝐼)] = ∫ ⋯ ⋯ ∫ 𝑑𝑥1 ⋯ 𝑑𝑥𝑁 max{𝐴(0, 𝑇) − 𝐾, 0} ∙
−∞
1
√(2𝜋𝜎 2 ∆𝑡)𝑁
−∞
𝑒
1
(− 2 ∑𝑁
(𝑥 −(𝑥𝑘−1 +𝜇 ∆𝑡))2
2𝜎 ∆𝑡 𝑘=1 𝑘
1
∆𝑡
𝜎 2 ∆𝑡(1 − 𝜏)𝜏
(𝑥𝑛 −𝑥𝑛−1 )𝜏+𝑥𝑛−1
𝐴(0, 𝑇) = ∫ 𝜔(𝜏)𝑒
(1 +
)𝑑𝜏
𝑇
2
0
Тогда алгоритм вычисления цены азиатского опциона методом Монте-Карло
с помощью континуального интеграла будет выглядеть следующим образом:
𝜎2
1.
Вычислить 𝑥0 = ln 𝑆0 , 𝜇 = 𝑟 −
2.
Сгенерировать значения 𝑥1 , 𝑥2 , … , 𝑥𝑛+1 со следующими
плотностями для 𝑥𝑖 :
16
2
.
1
√2𝜋𝜎 2 ∆𝑡)
𝑒
1
(− 2 (𝑥𝑖 −(𝑥𝑖−1 +𝜇∆𝑡))2
2𝜎 ∆𝑡
3. Имея значения 𝑥0 , 𝑥1 , … , 𝑥𝑛+1 , вычислить
𝑛+1 1
∆𝑡
𝜎 2 ∆𝑡(1 − 𝜏)𝜏
(𝑥𝑘 −𝑥𝑘−1 )𝜏+𝑥𝑘−1
𝐴(0, 𝑇) = ∑ ∫ 𝜔(𝜏)𝑒
(1 +
)𝑑𝜏
𝑇
2
𝑘=1 0
4. Повторить 1 – 3 шаги 𝑁 раз.
5. С помощью полученных на 4 шаге значений вычислить:
𝐶𝑗̅ = max(𝐴(0, 𝑇) − 𝐾, 0) , 𝑗 = 1, … , 𝑁
6. Найти значение цены опциона:
𝑁
𝑒 −𝑟𝑇
𝐶(𝑆0 , 𝑇) =
∑ 𝐶𝑗̅
𝑁
𝑖=1
7. Повторять все шаги до тех пор, пока значение, полученное на текущем
шаге, 𝐶(𝑆0 , 𝑇) не будет отличаться от значения, полученного на
предыдущем шаге, на достаточно малую величину.
17
Глава 3. Обзор вычислительных систем
Под вычислительной системой обычно понимают совокупность
аппаратного и программного обеспечения, которая предназначена для
решения задач пользователя. При этом система имеет вычислители, с
помощью которых можно осуществлять параллельную обработку.
В
сфере
суперкомпьютерных
технологий
задача
повышения
производительности систем всегда оставалась важной. В процессе развития
этих систем, принимались различные решения для достижения этой цели.
В 60-е годы появились интегральные микросхемы, что привело к
скачку
производительности
за
счет
принципов
конвейеризации
и
суперскалярности. На их основе были созданы микропроцессоры (CPU), а
позже были созданы графические ускорители (GPU). [9]
Большое
количество
параллельных
вычислений
и
сложность
масштабирования многопроцессорных и многоядерных систем привели к
использованию GPU не только для ускорения трехмерной графики, но и для
решения задач, которые обладают высокой степенью параллелизма. Позже
это привело к появлению вычислительной техники GPGPU. [14]
GPU состоит из однородных вычислительных элементов с общей
памятью, каждый из которых может исполнять тысячи потоков. Эти потоки
могут быть сгруппированы в блоки, которые имеют общий кэш и быструю
разделяемую память.
Так как современные GPU имеют высокую скорость доступа к модулям
памяти,
обработка
больших
массивов
данных
может
происходить
параллельно, при этом производительность достигает высоких значений.
Если изначально CPU и GPU были созданы для определенного класса
задач, а системы являлись гомогенными, то есть состояли из одного или
18
нескольких вычислителей одинаковой архитектуры, то сейчас часто
используют
гетерогенные
системы.
Гетерогенные
системы
могут
использовать универсальный процессор CPU и графический GPU совместно.
Стандартной гетерогенной системой является совокупность одного
CPU и одного или более GPU. Однако GPU используется как сопроцессор к
центральному процессору, который является «хостом», и называется
«устройством».
GPGPU (General Purpose computing for GPU) – это техника
использования GPU для расчетов, которые обычно выполняются на CPU. [10]
Существуют
различные
платформы
для
GPGPU.
Наиболее
популярными являются:
OpenCL
Nvidia CUDA
DirectCompute
AMD FireStream
И др.
3. 1. Обзор Nvidia CUDA
Компания Nvidia разработала программно-аппаратную архитектуру
параллельных вычислений CUDA (Compute Unified Device Architecture),
которая позволяет реализовать неграфические вычисления на графических
процессорах. [11][12][13][15]
CUDA-программа использует и CPU, и GPU. На CPU выполняется
последовательная часть кода. На GPU - параллельные участки кода,
выполняемые одновременно несколькими потоками.
Платформа CUDA предоставляет набор расширений для языков C/C++.
19
CUDA позволяет по собственному усмотрению разработчика управлять
доступом к набору инструкций и памяти.
С помощью CUDA можно определить так называемые ядра (kernel) –
функции, выполняющиеся на графическом процессоре. Ядро выполняется
такое же количество раз, сколько запускается потоков на GPU. То есть поток
выполняет ядро.
При этом ядро может быть вызвано из кода, выполняющегося на CPU,
но даже после вызова, код продолжит выполняться. Поэтому нужно
производить синхронизацию.
В соответствии с CUDA ядро выполняется на так называемом гриде,
который состоит из блоков, которые в свою очередь состоят из нитей
(потоков).
Разделение
на
блоки
необходимо
для
эффективного
взаимодействия нитей между собой в одном блоке, что дает некоторые
преимущества для параллельных алгоритмов.
CUDA имеет несколько видов памяти:
1.
Регистровая память. Регистры распределяются по нитям блока на
этапе компиляции и являются доступными им для чтения и записи на
все время исполнения ядра.
2.
Разделяемая память. Размещена на потоковом мультипроцессоре.
Распределяется по блокам, предоставляя каждому одно и то же
количество разделяемой памяти. Является самой быстрой, но в то же
время ограниченной.
3.
Локальная память. Размещена в DRAM (динамической памяти с
произвольным
доступом).
Используется
нитями
для
хранения
локальных данных.
4.
Глобальная память. Обычная память DRAM. Расположена на
GPU. Является самой медленной из всех типов памяти CUDA.
Выделяется на грид, то есть любая нить может пользоваться ей для
20
чтения и записи. Обычно используется для хранения данных большого
объема.
5.
Константная и текстурная память. Размещены в DRAM.
Обладают независимым кэшем. Имеют высокую скорость доступа.
Доступна всем нитям GPU для чтения. CPU может проводить как
запись, так и чтение.
21
Глава 4. Реализация с помощью CUDA
Для того чтобы использование GPU было эффективным, необходимо
использовать
некоторые
оптимизированные
под
него
параллельные
алгоритмы. Рассмотрим два алгоритма обработки массивов.
4. 1. Алгоритм параллельной редукции
Опр. 4. Редукцией массива 𝑥0 , 𝑥1 , 𝑥2 , … , 𝑥𝑛−1
называется следующее
выражение:
𝐴 = (((𝑥0 + 𝑥1 ) + 𝑥2 + ⋯ + 𝑥𝑛−1 )
Реализация последовательной редукции тривиальна.
Алгоритм параллельной редукции:
1.
Используется 𝑁 потоков.
2.
Каждым потоком вычисляются суммы двух элементов, которые
записываются в разделяемую память.
3.
Второй шаг повторяется 𝑙𝑜𝑔𝑁 раз, при этом уменьшая число
активных потоков в 2 раза на каждой итерации, также необходима
синхронизация после каждого шага.
На рисунке 1 на примере последовательности из 16 чисел наглядно
показан данный алгоритм. Подробное описание алгоритма можно найти в
данной книге [11].
4. 2. Алгоритм префиксных сумм
Рассмотрим алгоритм префиксных сумм.
22
Опр. 5. Префиксная сумма (prefix sum, scan) последовательности чисел
𝑥0 , 𝑥1 , 𝑥2 , …
-
это
последовательность
чисел
𝑦0 , 𝑦1 , 𝑦2 , …,
которые
вычисляются следующим образом:
𝑦0 = 𝑥0
𝑦1 = 𝑥0 + 𝑥1
𝑦2 = 𝑥0 + 𝑥1 + 𝑥2
…
Рис. 1.
Реализация последовательного алгоритма префиксных сумм (который
может запускаться на одном потоке на CPU, например) тривиальна.
Рассмотрим параллельную реализацию. Она стоит из двух частей, в
каждой их которых повторяется 𝑙𝑜𝑔2𝑁 шагов. На первом шаге строится
дерево, состоящее из сумм (так же, как в случае с параллельной редукцией),
на втором – результирующий массив по данному дереву сумм.
На рисунке 2 можно увидеть, как работает алгоритм для входной
последовательности из 16 чисел. Подробное описание алгоритма можно
увидеть в данной книге [11].
23
4. 3. Генерация случайных чисел.
Библиотека генерации случайных чисел cuRAND обеспечивает
высокопроизводительное GPU-ускоренное создание случайных чисел. [16]
Существуют четыре алгоритма генерации случайных чисел:
MRG32k3a
Генератор «Вихрь Мерсенна»
Генератор псевдослучайных чисел XORWOW
Генераторы квазислучайных чисел Соболя, включая поддержку
зашифрованных и 64-битных генераторов случайных чисел
Первые три являются псевдослучайными генераторами, последний –
квазислучайным.
Рис. 2.
24
При использовании методов Монте-Карло важную роль играют
качество и производительность генераторов случайных чисел.
Генератор псевдослучайных чисел создает все числа из диапазона с
равной вероятностью, поэтому генерация случайной величины не влияет на
возможность повторения ее генерации на следующем шаге. Генератор
квазислучайных чисел, в свою очередь, уменьшает вероятность ее
повторного появления, поэтому этот генератор покрывает пространство
равномернее.
В силу особенностей двух типов генераторов, представленных выше,
для
симуляции
Монте-Карло
выгоднее
квазислучайных чисел Соболя.
25
использовать
генераторы
Глава 5. Численные эксперименты.
Для запуска параллельной реализации алгоритма использованы три
GPGPU NVidia Tesla M2050. Этот графический процессор имеет следующие
характеристики:
Производительность операций с плавающей запятой
Производительность
операций
с
плавающей
515 ГФлоп
запятой 1.03 Тфлоп
одинарной точности (пиковая)
Полный объем специальной памяти
3 ГБ GDDR5
Макс. потребление энергии
225 Вт
Системный интерфейс
PCIe x16 Gen2
Количество ядер CUDA
448
В данной главе будет представлено сравнение точности и скорости
вычислений цены опциона на CPU и GPGPU. А также будет рассмотрена
возможность масштабирования.
Вычисления произведены для опциона типа «колл» со следующими
начальными данными:
1.
𝑆0 = 100
2.
𝐾 = 110
3.
𝑇 = 1 год
4.
𝜎 = 0.3
5.
𝑟 = 0.1
Если подставить данные значения в формулу Блэка-Шоулза (1)
получим следующее:
100
0.32
ln (
+ (0.1 +
)
110)
2 = 0.1656 ≈ 0.17
𝑑1 =
0.3
26
𝑑2 =
ln (
100
0.32
+ (0.1 −
)
)
110
2 = −0.1344 ≈ −0.13
0.3
𝑁(𝑑1 ) = 𝑁(0.17) = 0.5675
𝑁(𝑑2 ) = 𝑁(−0.13) = 1 − 𝑁(0.13) = 1 − 0.5517 = 0.4483
𝐶(𝑆0 , 0) = 100 ∗ 0.5675 − 110 ∗ 𝑒 −0.1 ∗ 0.4483 = 56.75 − 44.62 = 12.13
То есть цена такого европейского опциона будет равна 12,13. Стоит
отметить,
что
формула
Блэка-Шоулза
следует
определенным
предположениям, которые практически не могут соответствовать реальным
условиям.
5.1. Вычисление цены европейского опциона с помощью
стохастического дифференциального уравнения
Оценка европейского опциона методом Монте-Карло
18,0000
16,0000
14,0000
10,0000
Монте-Карло
8,0000
Блек-Шоулз
6,0000
4,0000
2,0000
0,0000
10
90
170
250
330
410
490
570
650
730
810
890
970
1050
1130
1210
1290
1370
1450
1530
1610
1690
1770
1850
1930
Цена опциона
12,0000
Количество случайных чисел
Рис. 3.
27
На графике рисунка 3 красной линией обозначено значение, которое
вычислено с помощью формулы Блэка-Шоулза. Было сгенерировано малое
количество случайных чисел, поэтому сходимость к значению, полученному
по формуле Блэка-Шоулза, недостаточная. Попробуем увеличить количество
случайных чисел.
На рисунке 4 видно, что при 20 млн. случайных величин сходимость
намного выше, чем в предыдущем эксперименте с точностью до двух знаков
после запятой.
Оценка европейского опциона методом Монте-Карло
12,2200
12,2000
12,1800
Цена опциона
12,1600
12,1400
12,1200
12,1000
Монте-Карло
12,0800
Блек-Шоулз
12,0600
12,0400
100 000
900 000
1 700 000
2 500 000
3 300 000
4 100 000
4 900 000
5 700 000
6 500 000
7 300 000
8 100 000
8 900 000
9 700 000
10 500 000
11 300 000
12 100 000
12 900 000
13 700 000
14 500 000
15 300 000
16 100 000
16 900 000
17 700 000
18 500 000
19 300 000
12,0200
Количество случайных чисел
Рис. 4.
Важно определить ускорение на GPU, и сравнить скорость вычислений
со скоростью на CPU.
На рисунке 5 показано ускорение на GPU, при различных наборах
случайных чисел. При увеличении числа случайных величин ускорение
растет, так как CPU исполняет поток инструкций последовательно, хоть и с
28
максимальной производительностью, а GPU исполняет большое число
потоков инструкций параллельно.
В таблице 1 представлен сравнительный анализ скорости вычислений
на CPU и GPU.
Кол-во случ. чисел
Время на CPU, мс
Время на GPU, мс
358 400
59
12
3 942 400
653
16
7 526 400
1242
20
11 110 400
1839
22
14 694 400
2442
26
20 428 800
3376
32
Табл. 1.
Таким образом, при 300 тыс. чисел ускорение равно 4.7, а при 2 млн. – 106.
Ускорение вычислений на GPU
120
100
60
CPU/GPU
40
20
0
358 400
1 075 200
1 792 000
2 508 800
3 225 600
3 942 400
4 659 200
5 376 000
6 092 800
6 809 600
7 526 400
8 243 200
8 960 000
9 676 800
10 393 600
11 110 400
11 827 200
12 544 000
13 260 800
13 977 600
14 694 400
15 411 200
16 128 000
16 844 800
17 561 600
18 278 400
18 995 200
19 712 000
20 428 800
Ускорение
80
Количество случайных чисел
Рис. 5.
29
5. 2. Вычисление цены азиатского опциона с помощью
стохастического дифференциального уравнения
Сложность азиатского опциона состоит в том, что он может быть
исполнен в любой момент до определенного момента, указанного в
контракте. Поэтому необходимо разделить время на отрезки и для каждой
точки считать цену базового актива и среднюю.
Используем те же начальные данные, которые указаны в начале главы.
На рисунке 6 видно среднее значение опциона, оно равно 3,92. Сходимость
цены к этому значению представлена на графике.
Рисунок 7 изображает то, как изменяется скорость вычислений цены
азиатского опциона в зависимости от использования CPU или GPU.
Можно заметить, что по сравнению с европейским опционом
ускорение намного ниже. Это связано с особенностями опционов азиатского
типа,
так
как
на
каждой
итерации
было
необходимо
вычислить
промежуточное значение цены базового актива и использовать алгоритм
префиксных сумм.
Так же, как и в предыдущем пункте, проведем сравнительный анализ
скорости вычислений на CPU и GPU.
Кол-во случ. чисел
Время на CPU, мс
Время на GPU, мс
358 400
56
142
3 942 400
616
327
7 526 400
1182
635
11 110 400
1749
923
14 694 400
2299
1204
20 428 800
3180
1696
Табл. 2.
В соответствии с таблицей 2, ускорение в среднем равно 2. Причины такого
малой разницы в скорости вычислений указаны выше.
30
358400
1075200
1792000
2508800
3225600
3942400
4659200
5376000
6092800
6809600
7526400
8243200
8960000
9676800
10393600
11110400
11827200
12544000
13260800
13977600
14694400
15411200
16128000
16844800
17561600
18278400
18995200
19712000
20428800
Время, мс
102400
716800
1331200
1945600
2560000
3174400
3788800
4403200
5017600
5632000
6246400
6860800
7475200
8089600
8704000
9318400
9932800
10547200
11161600
11776000
12390400
13004800
13619200
14233600
14848000
15462400
16076800
Цена опциона
Оценка азиатского опциона методом Монте-Карло
5
4,5
4
3,5
Цена опциона
Среднее значение
3
2,5
Количество случайных чисел
Рис. 6.
Скорость вычислений на GPU и CPU
3500
3000
2500
2000
1500
GPU
1000
CPU
500
0
Количество случайных величин
Рис. 7.
31
5. 3. Вычисление цены европейского опциона с помощью
континуального интеграла
Используем
те
же
начальные
данные
и
цену,
полученную
аналитически, равную 12,13.
Европейски опцион (континуальный интеграл)
14
13,5
Цена опциона
13
12,5
12
11,5
GPU
11
Среднее значение
10,5
10
358 400
2 508 800
4 659 200
6 809 600
8 960 000
11 110 400
13 260 800
15 411 200
17 561 600
19 712 000
21 862 400
24 012 800
26 163 200
28 313 600
30 464 000
32 614 400
34 764 800
36 915 200
39 065 600
41 216 000
43 366 400
45 516 800
47 667 200
49 817 600
9,5
Количество случайных величин
Рис. 8.
На рисунке 8 видно, что вычисляемая цена опциона сходится к
значению, вычисленному аналитически. Однако для достижения точности до
второго знака после запятой требуется большее количество случайных чисел,
чем в случае со стохастическим дифференциальным уравнением.
Сравним вычисления на CPU и GPU.
На рисунке 9 видно, что максимальное ускорение равно 74. Немного
ниже, чем в случае со стохастическим дифференциальным уравнением, но
это связано с некоторыми сложностями в данном алгоритме.
32
Ускорение вычислений на GPU
80
70
Ускорение
60
50
40
30
CPU/GPU
20
10
368 640
2 211 840
4 055 040
5 898 240
7 741 440
9 584 640
11 427 840
13 271 040
15 114 240
16 957 440
18 800 640
20 643 840
22 487 040
24 330 240
26 173 440
28 016 640
29 859 840
31 703 040
33 546 240
35 389 440
37 232 640
39 075 840
40 919 040
42 762 240
44 605 440
46 448 640
48 291 840
50 135 040
0
Количество случайных чисел
Рис. 9.
5. 4. Вычисление цены азиатского опциона с помощью
континуального интеграла
Используем те же начальные данные. Вычислим цену азиатского
опциона через континуальный интеграл.
На рисунке 10 видно, что вычисляемое значение опциона сходится к
значению 3,61. Разница с результатом работы алгоритма со стохастическим
уравнением (3,92) связана с основами подходов, однако является не
существенной.
33
Азиатский опцион (континуальный интеграл)
4,3
4,1
Цена опциона
3,9
3,7
3,5
Цена опциона
3,3
Среднее значение
3,1
29 859 840
28 385 280
26 910 720
25 436 160
22 487 040
23 961 600
21 012 480
19 537 920
18 063 360
16 588 800
15 114 240
13 639 680
12 165 120
9 216 000
10 690 560
7 741 440
6 266 880
4 792 320
3 317 760
368 640
1 843 200
2,9
Количество случайных чисел
Рис. 10.
На рисунке 11 можно увидеть, что максимальное ускорение равно 83.
Разница с ускорением вычислений с помощью стохастического уравнения
связана с тем, что не используется алгоритм префиксных сумм.
5. 5. Масштабируемость
Масштабируемость в параллельных вычислениях играет важную роль.
Так как при добавлении ресурсов важно чтобы производительность
действительно увеличивалась.
Проверим алгоритм с использованием континуального интеграла на
масштабируемость. Запустим его на одном, двух, трех GPU.
34
368 640
2 211 840
4 055 040
5 898 240
7 741 440
9 584 640
11 427 840
13 271 040
15 114 240
16 957 440
18 800 640
20 643 840
22 487 040
24 330 240
26 173 440
28 016 640
29 859 840
31 703 040
33 546 240
35 389 440
37 232 640
39 075 840
40 919 040
42 762 240
44 605 440
46 448 640
48 291 840
50 135 040
Ускорение
368 640
2 211 840
4 055 040
5 898 240
7 741 440
9 584 640
11 427 840
13 271 040
15 114 240
16 957 440
18 800 640
20 643 840
22 487 040
24 330 240
26 173 440
28 016 640
29 859 840
31 703 040
33 546 240
35 389 440
37 232 640
39 075 840
40 919 040
42 762 240
44 605 440
46 448 640
48 291 840
50 135 040
Ускорение
Ускорение вычислений на GPU
90
80
70
60
50
40
30
GPU
20
10
0
Количество случайных чисел
Рис. 11.
Маштабирование GPU
200
180
160
140
120
100
80
GPU
60
2 GPU
40
3 GPU
20
0
Количество случайных чисел
Рис. 12.
35
Рисунок 12 показывает хорошую масштабируемость данной задачи.
Помимо этого, каждое ускорение стремится к конкретному значению.
Например, при использовании одного GPU ускорение равно 72, при двух
GPU – 120, при трех GPU – 180.
Такие результаты связаны с тем, что сами вычисления крайне мало
зависят друг от друга.
36
Выводы
В данной работе проанализирована эффективность применения
технологии CUDA с использованием GPGPU NVIDIA для расчета цен
европейского и азиатского опционов методами Монте-Карло с помощью
стохастического дифференциального уравнения и континуального интеграла.
Для решения этой задачи с помощью симуляции Монте-Карло на GPU
использовались алгоритмы параллельной редукции и префиксных сумм, в
зависимости от типа алгоритма. Также был проведен анализ генераторов
случайных чисел и выбор оптимального из них для поставленной задачи.
Проведенные расчеты показали, что использование GPU достаточно
эффективно в применении к данной задаче. Для того чтобы достичь высокой
точности вычислений, требовалось сгенерировать большое количество
случайных величин, поэтому использование последовательного алгоритма на
СPU является менее эффективным. Реализация методов Монте-Карло на
GPGPU NVidia Tesla M2050 позволила достичь следующих результатов:
106-кратное ускорение для вычислений с помощью стохастического
дифференциального уравнения для европейских опционов
2-кратное ускорение для вычислений с помощью стохастического
дифференциального уравнения для азиатских опционов
70-кратное ускорение для вычислений с помощью континуального
интеграла для европейских опционов
75-кратное ускорение для вычислений с помощью континуального
интеграла для азиатских опционов
Вычисление цены азиатских опционов с помощью СДУ показало
неэффективность применения GPU, так как использовался алгоритм
префиксных
сумм.
Также
алгоритмы
масштабируемостью.
37
обладают
хорошей
Заключение
Так как, в настоящее время появляется больше суперкомпьютеров с
гетерогенной архитектурой, которые используют CPU и GPU одновременно,
технология CUDA применяется к огромному количеству задач, в том числе и
в финансовой математике.
Для вычисления цены опционов широко используются методы МонтеКарло. Эти алгоритмы хорошо подходят для реализации на GPU, так как они
основываются на большом количестве независимых операций, которые затем
помогают получить конечный результат. То есть методы Монте-Карло
обладают высокой степенью параллелизма.
По этим причинам была проведена данная работа. Сформулированные
задачи в рамках этой работы были выполнены, а ее результаты были
проанализированы.
38
Список литературы
[1]. F. Black, M. Scholes. The pricing of options and corporate liabilities // Journal
Political Economy, 1973. Vol. 81. No. 3. P. 637-659
[2]. Paul Glasserman, Monte Carlo Methods in Financial Engineering, Springer,
2003, 599 p.
[3]. Соболь И.М., Численные методы Монте-Карло. М.: Наука, 1973. 307 с.
[4]. Опцион, Википедия. https://ru.wikipedia.org/wiki/Опцион
[5]. Джон К. Халл. Опционы, фьючерсы и другие производные финансовые
инструменты. Изд. дом. Вильямс, 2008. 1044 c.
[6]. Hongbin Zhang. Pricing asian options using Monte Carlo methods //
Examensarbete i Matematik, U.U.D.M. Project Report. 2009 Vol. 2009:7
[7]. Vadim Linetsky. 14 approach to financial modeling and options pricing //
Computational Economics, 1998. Vol. 11: P. 129-163
[8]. Guido Montagna, Oreste Nicrosini. A path integral way to option pricing //
Physica A: Statistical Mechanics and its Applications, 2002. Volume 310, Issues 34, P. 450-466
[9]. Э. Таненбаум, М. ван Стеен. Распределенные системы. Принципы и
парадигмы. СПб.: Питер, 2003. 877 с.
[10]. Hyesoon Kim, Richard Vuduc, Sara Baghsorkhi. Performance Analysis and
Tuning for General Purpose Graphics Processing Units (GPGPU) // Morgan &
Claypool Publishers, 2012
[11]. А. В. Боресков и др. Предисл.: В. А. Садовничий. Параллельные
вычисления на GPU. Архитектура и программная модель CUDA: Учебное
пособие. Изд-во Московского университета, 2012, 336 стр.
[12]. А. В. Боресков, А. А. Харламов. Основы работы с технологией CUDA.
Изд. дом. ДМК Пресс, 2010, 232 стр.
[13]. Официальный сайт Nvidia
http://www.nvidia.ru/object/gpucomputingapplications-ru.html
39
[14]. Simoes B. General-purpose computing on the GPU (GPGPU).
http://www.think-techie.com/2009/09/general-purpose-computing-on-gpugpgpu.html
[15]. NVIDIA CUDA C Programming Guide. Version 4.2
http://developer.nvidia.com/nvidia-gpu-computing-documentation
[16].Библиотека cuRand http://www.nvidia.ru/object/tesla-gpu-acceleratedlibraries-curand-ru.html
40
Отзывы:
Авторизуйтесь, чтобы оставить отзыв