Python王牌加速库：深度学习下的障碍期权定价！

蒙特卡罗模拟需要数以百万计的路径来得到精确的答案，这需要大量的计算。Ryan等人得研究表明，可以训练深度学习模型对衍生品进行估值。深度学习模型是准确和快速的，能够产生比传统模型快一百万倍的估值。在今天的推文中，我们将使用一个 全连接网络来学习亚式障碍期权的定价模式 。采用蒙特卡罗模拟作为训练的定价依据。我们使用与上一篇文章相同的亚式障碍期权模型，参数如下：

• T：到期如（年）
• S：现货（美元）
• K：Strike（美元）
• sigma：波动率（per）
• r：无风险利率（per）
• mu：Drift Rate（per）
• B：Barrier（美元）

下面的内容主要包括两个主题：

• 使用蒙特卡罗定价动态数据集训练期权定价神的经网络模型。
• 使用蒙特卡罗定价静态数据集训练期权定价神经网络模型并进行推断。

2

批处理数据生成

数据集是深度学习训练的重要组成部分。我们将修改之前的单一亚式障碍期权定价代码来处理一批障碍期权定价。

加载库：

import cupy
import numpy as np
import math
import time
import torch
cupy.cuda.set_allocator(None)
from torch.utils.dlpack import from_dlpack

批量障碍期权定价模拟的CuPy版本如下：

cupy_batched_barrier_option = cupy.RawKernel(r'''
extern "C" __global__ void batched_barrier_option(float *d_s,const float T,const float * K,const float * B,const float * S0,const float * sigma,const float * mu,const float * r,const float * d_normals,const long N_STEPS,const long N_PATHS,const long N_BATCH)
{unsigned idx =  threadIdx.x + blockIdx.x * blockDim.x;unsigned stride = blockDim.x * gridDim.x;unsigned tid = threadIdx.x;const float tmp3 = sqrt(T/N_STEPS);for (unsigned i = idx; iK[batch_id] ? running_average-K[batch_id] : 0.f); d_s[i] = tmp2 * payoff;}
}''', 'batched_barrier_option')

注意，参数（K, B, S0, sigma, mu, r）以批处理长度的数组形式传入。输出数组是一个1-D 的二维数组。第一个维度用于 Batch，第二个维度用于 Path。。

通过输入两组选项参数进行测试：

N_PATHS = 2048000
N_STEPS = 365
N_BATCH = 2
T = 1.0K = cupy.array([110.0, 120.0], dtype=cupy.float32)
B = cupy.array([100.0, 90.0], dtype=cupy.float32)
S0 = cupy.array([120.0, 100.0], dtype=cupy.float32)
sigma = cupy.array([0.35, 0.2], dtype=cupy.float32)
mu = cupy.array([0.15, 0.1], dtype=cupy.float32)
r =cupy.array([0.05, 0.05], dtype=cupy.float32)

把这一切放进一个简单的函数来启动1GPU内核。每个Path的期权价格是相应路径terminal值的平均值。这可以很容易地通过Cupy函数平均值（axis=1）计算出来

def batch_run():number_of_threads = 256number_of_blocks = (N_PATHS * N_BATCH - 1) // number_of_threads + 1randoms_gpu = cupy.random.normal(0, 1, N_BATCH*N_PATHS * N_STEPS, dtype=cupy.float32)output = cupy.zeros(N_BATCH*N_PATHS, dtype=cupy.float32)cupy.cuda.stream.get_current_stream().synchronize()s = time.time()cupy_batched_barrier_option((number_of_blocks,), (number_of_threads,),(output, np.float32(T), K, B, S0, sigma, mu, r,randoms_gpu, N_STEPS, N_PATHS, N_BATCH))v = output.reshape(N_BATCH, N_PATHS).mean(axis=1)cupy.cuda.stream.get_current_stream().synchronize()e = time.time()print('time', e-s, 'v',v)
batch_run()

time 0.013919591903686523 v [21.22405    0.8480416]

这将为66ms中的这两组期权参数生成21.22和0.848的期权价格。

它的工作效率很高，因此我们将构造一个OptionDataSet类来包装上面的代码，以便我们可以在Pytorch中使用它。对于每个下一个元素，生成指定范围内的均匀分布随机期权参数，启动GPU内核计算期权价格，通过DLPack将CuPy数组转换为带有zero-copy的Pytorch张量。请注意我们是如何实现iterable Dataset接口的：

class OptionDataSet(torch.utils.data.IterableDataset):def __init__(self, max_len=10, number_path = 1000, batch=2, threads=256,seed=15):self.num = 0self.max_length = max_lenself.N_PATHS = number_pathself.N_STEPS = 365self.N_BATCH = batchself.T = np.float32(1.0)self.output = cupy.zeros(self.N_BATCH*self.N_PATHS, dtype=cupy.float32) self.number_of_blocks = (self.N_PATHS * self.N_BATCH - 1) // threads + 1self.number_of_threads = threadscupy.random.seed(seed)def __len__(self):return self.max_lengthdef __iter__(self):self.num = 0return selfdef __next__(self):if self.num > self.max_length:raise StopIterationX = cupy.random.rand(self.N_BATCH, 6, dtype=cupy.float32)X = X * cupy.array([200.0, 0.99, 200.0, 0.4, 0.2, 0.2], dtype=cupy.float32)X[:, 1] = X[:, 0] * X[:, 1]randoms = cupy.random.normal(0, 1, self.N_BATCH * self.N_PATHS * self.N_STEPS, dtype=cupy.float32)cupy_batched_barrier_option((self.number_of_blocks,), (self.number_of_threads,), (self.output, self.T, cupy.ascontiguousarray(X[:, 0]), cupy.ascontiguousarray(X[:, 1]), cupy.ascontiguousarray(X[:, 2]), cupy.ascontiguousarray(X[:, 3]), cupy.ascontiguousarray(X[:, 4]), cupy.ascontiguousarray(X[:, 5]), randoms, self.N_STEPS, self.N_PATHS, self.N_BATCH))Y = self.output.reshape(self.N_BATCH, self.N_PATHS).mean(axis=1)self.num += 1return (from_dlpack(X.toDlpack()), from_dlpack(Y.toDlpack()))

将所有与Pytorch数据集相关的内容都放到一个名为cupy_dataset.py的文件中：

%%writefile cupy_dataset.py 
import cupy
import numpy as np
import torch
from torch.utils.dlpack import from_dlpack
cupy.cuda.set_allocator(None)cupy_batched_barrier_option = cupy.RawKernel(r'''
extern "C" __global__ void batched_barrier_option(float *d_s,const float T,const float * K,const float * B,const float * S0,const float * sigma,const float * mu,const float * r,const float * d_normals,const long N_STEPS,const long N_PATHS,const long N_BATCH)
{unsigned idx =  threadIdx.x + blockIdx.x * blockDim.x;unsigned stride = blockDim.x * gridDim.x;unsigned tid = threadIdx.x;const float tmp3 = sqrt(T/N_STEPS);for (unsigned i = idx; iK[batch_id] ? running_average-K[batch_id] : 0.f); d_s[i] = tmp2 * payoff;}
}''', 'batched_barrier_option')class OptionDataSet(torch.utils.data.IterableDataset):def __init__(self, max_len=10, number_path = 1000, batch=2, threads=256,seed=15):self.num = 0self.max_length = max_lenself.N_PATHS = number_pathself.N_STEPS = 365self.N_BATCH = batchself.T = np.float32(1.0)self.output = cupy.zeros(self.N_BATCH*self.N_PATHS, dtype=cupy.float32) self.number_of_blocks = (self.N_PATHS * self.N_BATCH - 1) // threads + 1self.number_of_threads = threadscupy.random.seed(seed)def __len__(self):return self.max_lengthdef __iter__(self):self.num = 0return selfdef __next__(self):if self.num > self.max_length:raise StopIterationX = cupy.random.rand(self.N_BATCH, 6, dtype=cupy.float32)X = X * cupy.array([200.0, 0.99, 200.0, 0.4, 0.2, 0.2], dtype=cupy.float32)X[:, 1] = X[:, 0] * X[:, 1]randoms = cupy.random.normal(0, 1, self.N_BATCH * self.N_PATHS * self.N_STEPS, dtype=cupy.float32)cupy_batched_barrier_option((self.number_of_blocks,), (self.number_of_threads,), (self.output, self.T, cupy.ascontiguousarray(X[:, 0]), cupy.ascontiguousarray(X[:, 1]), cupy.ascontiguousarray(X[:, 2]), cupy.ascontiguousarray(X[:, 3]), cupy.ascontiguousarray(X[:, 4]), cupy.ascontiguousarray(X[:, 5]), randoms, self.N_STEPS, self.N_PATHS, self.N_BATCH))Y = self.output.reshape(self.N_BATCH, self.N_PATHS).mean(axis=1)self.num += 1return (from_dlpack(X.toDlpack()), from_dlpack(Y.toDlpack()))

覆盖cupy_dataset.py

这里是一个测试代码样本，有10个数据点、 batch为 16：

ds = OptionDataSet(10, number_path=100000, batch=16, seed=15)
for i in ds:print(i[1])

我们可以实现相同的代码使用Numba加速计算在GPU：

import numba
from numba import cuda@cuda.jit
def batch_barrier_option(d_s, T, K, B, S0, sigma, mu, r, d_normals, N_STEPS, N_PATHS, N_BATCH):ii = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.xstride = cuda.gridDim.x * cuda.blockDim.xtmp3 = math.sqrt(T/N_STEPS)for i in range(ii, N_PATHS * N_BATCH, stride):batch_id = i // N_PATHSpath_id = i % N_PATHStmp1 = mu[batch_id]*T/N_STEPStmp2 = math.exp(-r[batch_id]*T)running_average = 0.0s_curr = S0[batch_id]for n in range(N_STEPS):s_curr += tmp1 * s_curr + sigma[batch_id]*s_curr*tmp3*d_normals[path_id + batch_id * N_PATHS + n * N_PATHS * N_BATCH]running_average = running_average + 1.0/(n + 1.0) * (s_curr - running_average)if i==0 and batch_id == 2:print(s_curr)if running_average <= B[batch_id]:breakpayoff = running_average - K[batch_id] if running_average > K[batch_id] else 0d_s[i] = tmp2 * payoffclass NumbaOptionDataSet(object):def __init__(self, max_len=10, number_path = 1000, batch=2, threads=512, seed=15):self.num = 0self.max_length = max_lenself.N_PATHS = number_pathself.N_STEPS = 365self.N_BATCH = batchself.T = np.float32(1.0)self.output = cupy.zeros(self.N_BATCH*self.N_PATHS, dtype=cupy.float32) self.number_of_blocks = (self.N_PATHS * self.N_BATCH - 1) // threads + 1self.number_of_threads = threadscupy.random.seed(seed)def __len__(self):return self.max_lengthdef __iter__(self):self.num = 0return selfdef __next__(self):if self.num > self.max_length:raise StopIterationX = cupy.random.rand(self.N_BATCH, 6, dtype=cupy.float32)X = X * cupy.array([200.0, 0.99, 200.0, 0.4, 0.2, 0.2], dtype=cupy.float32)X[:, 1] = X[:, 0] * X[:, 1]randoms = cupy.random.normal(0, 1, self.N_BATCH * self.N_PATHS * self.N_STEPS, dtype=cupy.float32)batch_barrier_option[(self.number_of_blocks,), (self.number_of_threads,)](self.output, self.T, X[:, 0], X[:, 1], X[:, 2], X[:, 3], X[:, 4], X[:, 5], randoms, self.N_STEPS, self.N_PATHS, self.N_BATCH)o = self.output.reshape(self.N_BATCH, self.N_PATHS)Y = o.mean(axis = 1) self.num += 1return (from_dlpack(X.toDlpack()), from_dlpack(Y.toDlpack()))
ds = NumbaOptionDataSet(10, number_path=100000, batch=16, seed=15)
for i in ds:print(i[1])

3

模型

为了将期权参数映射到价格，我们使用了6层全连接神经网络，其隐含维度为512。将此DL价格模型写入model.py： 

%%writefile model.py
import torch.nn as nn
import torch.nn.functional as F
import torchclass Net(nn.Module):def __init__(self, hidden=1024):super(Net, self).__init__()self.fc1 = nn.Linear(6, hidden)self.fc2 = nn.Linear(hidden, hidden)self.fc3 = nn.Linear(hidden, hidden)self.fc4 = nn.Linear(hidden, hidden)self.fc5 = nn.Linear(hidden, hidden)self.fc6 = nn.Linear(hidden, 1)self.register_buffer('norm',torch.tensor([200.0,198.0,200.0,0.4,0.2,0.2]))def forward(self, x):x = x / self.normx = F.elu(self.fc1(x))x = F.elu(self.fc2(x))x = F.elu(self.fc3(x))x = F.elu(self.fc4(x))x = F.elu(self.fc5(x))return self.fc6(x)

覆盖model.py

输入参数首先通过除以（200.0,198.0,200.0,0.4,0.2,0.2）缩小到0-1范围。然后在ELu激活函数后，将其5次隐射到隐藏维度512。选择ELu是因为我们需要计算参数的二阶微分。如果使用ReLu，二阶微分总是0。最后一层是线性层，它将隐藏维度映射到预测的期权价格。

在训练方面，我们使用了一个高级库Ignite来训练PyTorch中的神经网络：

https://github.com/pytorch/ignite

我们使用MSELoss作为损失函数，Adam作为优化器，CosineAnnealingScheduler作为学习率调度器。下面的代码将随机期权数据提供给定价模型进行训练。

from ignite.engine import Engine, Events
from ignite.handlers import Timer
from torch.nn import MSELoss
from torch.optim import Adam
from ignite.contrib.handlers.param_scheduler import CosineAnnealingScheduler
from ignite.handlers import ModelCheckpoint
from model import Net
from cupy_dataset import OptionDataSet
timer = Timer(average=True)
model = Net().cuda()
loss_fn = MSELoss()
optimizer = Adam(model.parameters(), lr=1e-3)
dataset = OptionDataSet(max_len=10000, number_path = 1024, batch=4800)def train_update(engine, batch):model.train()optimizer.zero_grad()x = batch[0]y = batch[1]y_pred = model(x)loss = loss_fn(y_pred[:,0], y)loss.backward()optimizer.step()return loss.item()trainer = Engine(train_update)
log_interval = 100scheduler = CosineAnnealingScheduler(optimizer, 'lr', 1e-4, 1e-6, len(dataset))
trainer.add_event_handler(Events.ITERATION_STA

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！