Python王牌加速库：奇异期权定价的利器！

1

前言

在金融领域，计算效率有时可以直接转化为交易利润。 量化分析师面临着在研究效率和计算效率之间进行权衡的挑战 。使用Python可以生成简洁的研究代码，从而提高了研究效率。但是，一般的Python代码速度很慢，不适合用于生产环境。在这篇文章中，我们将探索如何使用 Python的GPU库来高性能实现奇异期权定价领域遇到的问题 。

2

定价计算概述

Black-Scholes模型可以有效地用欧洲行权规则为plain vanilla定价。像障碍（Barrier）期权和篮子（Basket ）期权这样的期权具有复杂的结构。蒙特卡罗模拟是一种有效的定价方法。为了得到一个精确的价格和一个小的变动，你需要许多模拟路径，计算十分密集。

幸运的是，每个模拟路径都是独立的， 大家可以利用多核NVIDIA GPU在一个节点内加速计算，甚至在必要时将其扩展到多个服务器 。由于独立路径的并行化，使用GPU可以将计算速度提高几个数量级。

传统上，对GPU的蒙特卡罗仿真是在CUDA C/ C++代码中实现的。大家必须明确地管理内存并编写大量样板代码，这对代码维护和生产效率提出了挑战。

最近，Deep Learning Derivatives（Ryan et al，2018）的论文被引入到使用深度神经网络来近似期权定价模型。该方法利用计算时间与推理时间进行定价训练，与GPU上的蒙特卡罗模拟相比，它实现了额外的数量级加速，这使得在生产环境中的实时奇异期权定价成为一个现实目标。

在这篇文章中介绍的方法对奇异期权类型没有任何限制。它适用于任何可以用蒙特卡罗方法模拟的期权定价模型。

在不失一般性的情况下，大家可以使用 亚式障碍期权 作为一个示例。亚式障碍期权是亚式期权和障碍期权的混合。衍生品价格取决于标的资产价格S、执行价格K和障碍价格B的平均值。以上下看涨期权离散化亚洲障碍期权为例。

• 如果标的资产的平均价格低于这一水平，则该期权无效。

• 资产现货价格S通常在建模中被认为是属于几何布朗运动，它有三个参数：现货价格、波动率和漂移率。

• 期权的价格是到期时的预期利润相对于当前价值的折现。

• 期权的路径依赖性使得对期权价格的解析解成为不可能。

这是使用蒙特卡罗模拟定价的一个很好的示例。你需要一个至少16GB的GPU来复现这个结果。

3

第1部分：使用GPU Python库进行蒙特卡洛定价

NVIDIA GPU被设计用来使用大量线程进行并行计算。蒙特卡罗仿真是在GPU中可以很好加速的算法之一。在下面的小节中，大家将看到在传统的CUDA代码中使用蒙特卡罗模拟，然后在Python中使用不同的库实现相同的算法。

CUDA方法

传统上，蒙特卡罗期权定价是在CUDA C/ C++中实现的。下面的CUDA C/ C++代码示例使用蒙特卡罗方法计算期权价格：

#include <vector>
#include <stdio.h>
#include <iostream>
#include <chrono>
#include <cuda_runtime.h>
#include <helper_cuda.h>
#include <curand.h>

#define CHECKCURAND(expression)                         \
  {                                                     \
    curandStatus_t status = (expression);                         \
    if (status != CURAND_STATUS_SUCCESS) {                        \
      std::cerr << "Curand Error on line " << __LINE__<< std::endl;     \
      std::exit(EXIT_FAILURE);                                          \
    }                                                                   \
  }

// atomicAdd is introduced for compute capability >=6.0
#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600
#else
__device__ double atomicAdd(double* address, double val)
{
      printf("device arch <=600\n");
        unsigned long long int* address_as_ull = (unsigned long long int*)address;
          unsigned long long int old = *address_as_ull, assumed;
            do {
                    assumed = old;
                        old = atomicCAS(address_as_ull, assumed,
                                                    __double_as_longlong(val + __longlong_as_double(assumed)));
                          } while (assumed != old);
              return __longlong_as_double(old);
}
#endif

__global__ void sumPayoffKernel(float *d_s, const unsigned N_PATHS, double *mysum)
{
  unsigned idx =  threadIdx.x + blockIdx.x * blockDim.x;
  unsigned stride = blockDim.x * gridDim.x;
  unsigned tid = threadIdx.x;

  extern __shared__ double smdata[];
  smdata[tid] = 0.0;

  for (unsigned i = idx; i<N_PATHS; i+=stride)
  {
    smdata[tid] += (double) d_s[i];
  }

  for (unsigned s=blockDim.x/2; s>0; s>>=1)
  {
    __syncthreads();
    if (tid < s) smdata[tid] += smdata[tid + s];
  }

  if (tid == 0)
  {
    atomicAdd(mysum, smdata[0]);
  }
}

__global__ void 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)
{
  unsigned idx =  threadIdx.x + blockIdx.x * blockDim.x;
  unsigned stride = blockDim.x * gridDim.x;
  const float tmp1 = mu*T/N_STEPS;
  const float tmp2 = exp(-r*T);
  const float tmp3 = sqrt(T/N_STEPS);
  double running_average = 0.0;

  for (unsigned i = idx; i<N_PATHS; i+=stride)
  {
    float s_curr = S0;
    for(unsigned n = 0; n < N_STEPS; n++){
       s_curr += tmp1 * s_curr + sigma*s_curr*tmp3*d_normals[i + n * N_PATHS];
       running_average += (s_curr - running_average) / (n + 1.0) ;
       if (running_average <= B){
           break;
       }
    }

    float payoff = (running_average>K ? running_average-K : 0.f);
    d_s[i] = tmp2 * payoff;
  }
}

int main(int argc, char *argv[]) {
  try {
    // declare variables and constants
    size_t N_PATHS = 8192000;
    size_t N_STEPS = 365;
    if (argc >= 2)  N_PATHS = atoi(argv[1]);

    if (argc >= 3)  N_STEPS = atoi(argv[2]);

    const float T = 1.0f;
    const float K = 110.0f;
    const float B = 100.0f;