深度学习在期权估值上的尝试

GhatGPT技术的推出，使得普通人可以用机器学习方法来解决问题：不需要精通编程技巧，甚至不需要深刻了解模型，这在以前几乎是不可以想象的(很多人开始焦虑了)。我有一些在衍生品估值方面的经验，最近开始接触深度学习，我觉得两个领域有些东西是相通的。比如，看涨期权的支付结构就是机器学习中的线性整流函(ReLU)，用多个看涨期权的组合来逼近任意的欧式支付结构就是机器学习中的通用近似定理，实现期权参数敏感度的快速计算的AAD(Adjoint Algorithmic Differentiation)算法无非就是深度学习中普遍用到的反向传播算法。

将机器学习应用在期权估值上，我想是很靠谱的，期权的估值与希腊字母在参数空间里是若干 (确定性的) 超曲面，根据通用近似定理，我们就可以用神经网络来实现映射？一旦实现，我们就可以摆脱对蒙特卡洛计算或有限差分算法(特别是高维)的漫长等待，取而代之的是经过AI芯片加速的矩阵乘法。

在下文中，我将把简单神经网络套用在欧式期权的估值上，分别针对Black-Scholes模型和局部波动率模型。这种做法在实际应用中可能没什么意义。这里一个原因是受限于计算资源，另一个原因，如标题所述，本文主要是一种尝试，希望能起到抛砖引玉的作用，我相信一些比较标准化的奇异期权（如障碍期权，亚式，经典雪球期权等可以按类似的思路去做，区别体现在输入参数的数目和神经网络的复杂度上)。

本文利用了现成的深度学习框架PyTorch，原始代码结构由大模型Kimi完成(要是闲得蛋疼，我可以用鲁迅或余华的文风给这篇文章润色)，训练目标值的准备是用QuantLib完成。比较偷懒，没什么图片，折腾出来的代码全部提供。废话不多说了，开始吧。

1. Black-Scholes模型

在Black-Scholes模型中，标的波动率$\sigma$为常数。

1.1 模型

对于欧式期权来说解析解就几乎可以解决所有问题，但在这里我们用它为例主要是为了展示模型(主要还是生成训练样本数据快)。一个欧式看涨期权，假设标的价格为$\text{S}$，行权价为$\text{K}$，无风险利率为$\text{r}$，年化到期期限为$\text{T}$，波动率为$\sigma$，则它的估值和几个常用的希腊字母计算如下

$V=S\cdot N(d_1)-K\cdot e^{-rT}\cdot N(d_2),$

$\Delta =N(d_1),$

$\Gamma =\frac{N^{\prime }(d_1)}{S\cdot \sigma \sqrt{T}},$

$\upsilon =S\cdot \sqrt{T}\cdot N^{\prime }(d_1),$

$\Theta =-\frac{S\cdot N^{\prime }(d_1)\cdot \sigma }{2\sqrt{T}}-K\cdot e^{-rT}\cdot N(d_2),$

$\rho =T\cdot K\cdot e^{-rT}\cdot N(d_2).$

$d_1$和$d_2$自行脑补。我们采用一个简单神经网络模型。为了减少输入参数，我们固定$\text{S}=1$，则行权价为相对值，最终的输入向量为$[K,r,T,\sigma ]$，输出向量为 $[V,\Delta ,\Gamma ,\upsilon ,\Theta ,\rho ]$。结合常用的参数范围，我们生成均匀的空间网格$linspace(0.5,1.5,101)\times linspace(0.0,0.1,11)\times linspace(0.01,2.01,201)\times linspace(0.01,1.01,101)$，这里$linspace$表示Numpy中的同名函数。以网格点为样本输入，筛选掉行权价在3个扩散标准差之外的样本，即样本参数满足限制

$(r-0.5{\sigma }^2)T-3\sigma \sqrt{T}<ln(K)<(r-0.5{\sigma }^2)T+3\sigma \sqrt{T}$，

剩下样本数约2005万，通过Black-Scholes公式计算出对应的输出。随机选择其中的90%为训练集，余下的为测试集。由于这是一个回归问题，我们选用带权重的MSE(mean square error)作为损失函数，即

$L=\frac{1}{N}{\sum }_{i=1}^N[w_V(V_i-{\hat{V}}_i{)}^2+w_{\Delta }({\Delta }_i-{\hat{\Delta }}_i{)}^2+w_{\Gamma }({\Gamma }_i-{\hat{\Gamma }}_i{)}^2+w_{\upsilon }({\upsilon }_i-{\hat{\upsilon }}_i{)}^2+w_{\Theta }({\Theta }_i-{\hat{\Theta }}_i{)}^2+w_{\rho }({\rho }_i-{\hat{\rho }}_i{)}^2]$

我们设置$[w_V,w_{\Delta },w_{\Gamma },w_{\upsilon },w_{\Theta },w_{\rho }]=[1.0,1.0,0.01,1.0,0.01,1.0]$(归一化前的值)，这样做的原因一方面是根据实际应用中关注度，另一方面是因为$\Gamma$和$\Theta$出现极端值的概率较高，我们希望减小其对训练效率的影响。

1.2 结果

选择的神经网络大小为$4\times 50\times 50\times 6$，这表示有两个50个神经元的隐藏层，学习率设置为0.001，在训练集上前后传播10次。最终，我们在测试集上的误差大约为$7.92\times 10^{-4}$，这是一个加权的结果，拆分到每一项的结果如下

-	$Weighted$	$V$	$\Delta$	$\Gamma$	$\upsilon$	$\Theta$	$\rho$
$MSE$	$7.92\times 10^{-4}$	$5.07\times 10^{-5}$	$1.58\times 10^{-4}$	$2.67\times 10^{-1}$	$1.23\times 10^{-4}$	$1.24\times 10^{-3}$	$1.74\times 10^{-4}$

不出意外，我们在权重较大的值上的训练效果较好。对于一般的期权，$5.07\times 10^{-5}$ 的估值方差并不是一个很精确的结果(对于前台报价来说)。如果进一步提高网络的复杂度，误差可能会更小，但是要注意避免过拟合。

1.3 改进

不需要那么痛苦。我们可以利用反向传播算法直接计算出$V$对于$S$，$r$， $T$， $\sigma$的导数，在这种情况下，模型的训练目标只有一个，就是$V$。直觉上，对于相同的网络深度和宽度，这样导致的训练效率会更高。
由于输入向量中不含有$S$，我们可能将输入向量扩充为$[1.0,K,r,T,\sigma ]$，但这是不可能的！考虑$S$从1变化到$1+ϵ$，如果重新归一化，输入向量中唯一会变化的就是$K$，变成$K/(1+ϵ)$，$\Delta$的计算就变成

$\Delta ={\lim }_{ϵ\to 0}\frac{(1+ϵ)\hat{V}(K/(1+ϵ))-\hat{V}(K)}{ϵ}=-K\cdot \frac{\partial \hat{V}(K)}{\partial K}+\hat{V}(K)$

$\Gamma$的计算需要涉及二阶导数，我们暂时不考虑。

PyTorch的神经网络计算框架中存储了模型参数，以及损失函数对于模型参数的导数，而我们需要的是输出值对于输入参数的导数，可以进行简单的运算转换得到。假设输入向量为$\mathbf{X}$，输出值为$y$，从输入层向第一个隐层有仿射变换

$\mathbf{Z}={\mathbf{W}}^T\mathbf{X}+b$

考虑单个样本的损失函数

$L(\hat{y},y)=\hat{y}-y$

两边对$\mathbf{X}$进行求导并参考链式法

$\frac{\partial L}{\partial \mathbf{X}}=\frac{\partial \hat{y}}{\partial \mathbf{X}},$

$\frac{\partial L}{\partial \mathbf{Z}}\mathbf{W}=\frac{\partial \hat{y}}{\partial \mathbf{X}}.$

左边需要的导数可以通过一次反向传播得到(不要更新模型参数)，和转移矩阵相乘就得到了需要的导数$\frac{\partial \hat{y}}{\partial \mathbf{X}}$。

$\Gamma$怎么办？难搞。回到一个基础的问题，我们是否要将希腊字母放入训练目标中？可能根本就不需要。第一，训练希腊字母会增加计算量，无论是参数数目上，还是收敛难度上(许多希腊字母呈现出各种陡峭的变化趋势)；第二，借助反向传播算法，我们可以快捷地计算出一些希腊字母(可以拿来吹牛逼的AAD)，如前所示；第三，希腊字母的定义多变，比如说隐含波动率曲面存在下的$\upsilon$，不同Sticky Rule下的$\Delta$，一天的$\Theta$，高阶值，等等。但是，只要我们有不同情景下的$V$(若干次神经网络的前向传播)，就可以通过“万金油”的差分办法算出它们，如果$V$算得足够精确，希腊字母值也可以算到不错的精度。因此，只关注对$V$的训练是一个不错的选择，在下文中，我们的目标值都只有$V$。

我们看一下数值结果，如下

-	$V$	$\Delta$	$\Gamma$	$\upsilon$	$\Theta$	$\rho$
$MS{E}_{full}$	$5.07\times 10^{-5}$	$1.58\times 10^{-4}$	$2.67\times 10^{-1}$	$1.23\times 10^{-4}$	$1.24\times 10^{-3}$	$1.74\times 10^{-4}$
$MS{E}_{aad}$	$7.66\times 10^{-7}$	$4.44\times 10^{-4}$	-	$3.98\times 10^{-4}$	$4.02\times 10^{-4}$	$9.33\times 10^{-4}$

不出意外，在相同的训练集和测试集下，只训练$V$得到的$MS{E}_{aad}$小了两个量级，通过AAD计算得到的希腊字母也能和直接训练保持相当的精度。抱歉这里没有用有限差分计算$\Gamma$，有兴趣的可以自行计算。

1.4 代码

代码如下，C++计算导出的CSV文件比较大就不提供了，很容易自己生成。我是Mac下运行的，自行切换CUDA。比较耗时间的训练，更耗时间的是逐个样本的通过AAD计算测试集，但也能控制在5个小时内，自行把握吧。

european_call_full.py

import math
import time

import numpy as np
from sklearn.model_selection import train_test_split
from torch import nn, optim
import torch
from torch.utils.data import TensorDataset, DataLoader
import scipy.stats as stats
import pandas as pd

# 设置设备
device = torch.device('mps' if torch.backends.mps.is_available() else 'cpu')


# Black-Scholes 欧式看涨期权模型
def black_scholes_call(K, r, T, sigma):
    S = 1  # 标的价格固定为1
    d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)

    call_price = S * stats.norm.cdf(d1) - K * math.exp(-r * T) * stats.norm.cdf(d2)
    delta = stats.norm.cdf(d1)
    gamma = stats.norm.pdf(d1) / (S * sigma * math.sqrt(T))
    vega = S * stats.norm.pdf(d1) * math.sqrt(T) / math.sqrt(2 * math.pi)
    theta = -S * gamma * sigma / (2 * math.sqrt(T)) - r * K * math.exp(-r * T) * stats.norm.cdf(d2)
    rho = K * math.exp(-r * T) * T * stats.norm.cdf(d2)

    return call_price, delta, gamma, vega, theta, rho


class WeightedMSELoss(nn.Module):
    def __init__(self, weight=None, reduction='mean'):
        super(WeightedMSELoss, self).__init__()
        self.weight = weight
        self.reduction = reduction

    def forward(self, input, target):
        # 计算MSE损失
        loss = (input - target) ** 2

        # 如果提供了权重，则应用权重
        if self.weight is not None:
            assert len(self.weight) == loss.size(-1), "Weight dimension must match the last dimension of the input"
            # loss = loss * self.weight
            loss = loss * self.weight

        # 根据reduction参数处理损失
        if self.reduction == 'none':
            return loss
        elif self.reduction == 'mean':
            return loss.mean()
        elif self.reduction == 'sum':
            return loss.sum()
        else:
            raise ValueError("Invalid reduction mode: {}".format(self.reduction))


# 构建神经网络模型
class OptionPricingNN(nn.Module):
    def __init__(self):
        super(OptionPricingNN, self).__init__()
        self.layers = nn.Sequential(
            nn.Linear(4, 50),
            nn.ReLU(),
            nn.Linear(50, 50),
            nn.ReLU(),
            nn.Linear(50, 6)
        )

    def forward(self, x):
        return self.layers(x)


if __name__ == "__main__":
    # K_range = np.linspace(0.5, 1.5, 101)
    # r_range = np.linspace(0.0, 0.1, 11)
    # T_range = np.linspace(0.01, 2.01, 201)
    # sigma_range = np.linspace(0.01, 1.01, 101)
    # 生成训练样本
    # K_grid, r_grid, T_grid, sigma_grid = np.meshgrid(K_range, r_range, T_range, sigma_range)
    # params = np.vstack((K_grid.ravel(), r_grid.ravel(), sigma_grid.ravel(), T_grid.ravel())).T
    # results = np.array([black_scholes_call(*param) for param in params])

    params = pd.read_csv('european/strike_rate_term_sigma.csv', index_col=0).values
    results = pd.read_csv('european/npv_delta_gamma_vega_theta_rho.csv', index_col=0).values

    filtered_rows = ((np.log(params[:, 0]) < (params[:, 1] - 0.5 * params[:, 3]**2) * params[:, 2]
                      + 3 * params[:, 3] * np.sqrt(params[:, 2]))
                   & (np.log(params[:, 0]) > (params[:, 1] - 0.5 * params[:, 3]**2) * params[:, 2]
                      - 3 * params[:, 3] * np.sqrt(params[:, 2])))
    params = params[filtered_rows]
    results = results[filtered_rows]
    print(results.shape)

    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(params, results, test_size=0.1, random_state=2024)

    # 转换数据为张量
    X_train = torch.tensor(X_train.astype(np.float32)).to(device)
    y_train = torch.tensor(y_train.astype(np.float32)).to(device)
    X_test = torch.tensor(X_test.astype(np.float32)).to(device)
    y_test = torch.tensor(y_test.astype(np.float32)).to(device)

    # 构建TensorDataset和DataLoader
    train_dataset = TensorDataset(X_train, y_train)
    test_dataset = TensorDataset(X_test, y_test)

    batch_size = 64
    train_loader = DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True)
    test_loader = DataLoader(dataset=test_dataset, batch_size=batch_size, shuffle=False)

    # 初始化模型、损失函数和优化器
    model = OptionPricingNN().to(device)
    weight = torch.tensor([1.0, 1.0, 0.01, 1.0, 0.01, 1.0]).to(device)
    weight = weight / weight.sum()
    criterion = WeightedMSELoss(weight=weight)
    optimizer = optim.Adam(model.parameters(), lr=0.001)

    print("Training begins: ")

    # 训练模型
    epochs = 10
    for epoch in range(epochs):
        model.train()  # 设置模型为训练模式
        running_loss = 0.0
        for batch_idx, (data, target) in enumerate(train_loader):
            optimizer.zero_grad()
            output = model(data)
            loss = criterion(output, target)
            loss.backward()
            optimizer.step()
            running_loss += loss.item()

        print(f'Epoch {epoch + 1}/{epochs}, Loss: {running_loss / len(train_loader)}')

    torch.save(model.state_dict(), 'european/model_weights_4_50_50_6.pth')

    #model.load_state_dict(torch.load('european/model_weights_4_50_50_6.pth', map_location=torch.device(device)))

    # 在测试集上评估模型
    print("Evaluation begins: ")
    start_time = time.time()  # 记录开始时间

    model.eval()  # 设置模型为评估模式
    criterion2 = WeightedMSELoss(weight=weight, reduction='sum')
    criterion3 = nn.MSELoss(size_average=None, reduce=None, reduction='sum')

    test_loss, npv_loss, delta_loss, gamma_loss, vega_loss, theta_loss, rho_loss = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

    with torch.no_grad():
        for data, target in test_loader:
            output = model(data)
            test_loss += criterion2(output, target).item()

            npv_loss += criterion3(output[:, [0]], target[:, [0]]).item()
            delta_loss += criterion3(output[:, [1]], target[:, [1]]).item()
            gamma_loss += criterion3(output[:, [2]], target[:, [2]]).item()
            vega_loss += criterion3(output[:, [3]], target[:, [3]]).item()
            theta_loss += criterion3(output[:, [4]], target[:, [4]]).item()
            rho_loss += criterion3(output[:, [5]], target[:, [5]]).item()

        print(f'Test  Loss: {test_loss / len(test_dataset)}')
        print(f'NPV   Loss: {npv_loss / len(test_dataset)}')
        print(f'Delta Loss: {delta_loss / len(test_dataset)}')
        print(f'Gamma Loss: {gamma_loss / len(test_dataset)}')
        print(f'Vega  Loss: {vega_loss / len(test_dataset)}')
        print(f'Theta Loss: {theta_loss / len(test_dataset)}')
        print(f'Rho   Loss: {rho_loss / len(test_dataset)}')

    end_time = time.time()  # 记录结束时间

    print("执行时间：", end_time - start_time)

european_call_aad.py

import time

import numpy as np
from sklearn.model_selection import train_test_split
from torch import nn, optim
import torch
from torch.utils.data import TensorDataset, DataLoader
import pandas as pd

# 设置设备
device = torch.device('mps' if torch.backends.mps.is_available() else 'cpu')


# 构建神经网络模型
class OptionPricingNN(nn.Module):
    def __init__(self, width=50):
        super(OptionPricingNN, self).__init__()
        self.width = width
        self.fc1 = nn.Linear(4, width)
        self.layers = nn.Sequential(
            nn.ReLU(),
            nn.Linear(width, width),
            nn.ReLU(),
            nn.Linear(width, 1)
        )

    def forward(self, x):
        x = self.fc1(x)
        x = self.layers(x)

        return x


class DiffLoss(nn.Module):
    def __init__(self):
        super(DiffLoss, self).__init__()

    def forward(self, input, target):
        # 计算MSE损失
        loss = input - target

        return loss


if __name__ == "__main__":
    params = pd.read_csv('european/strike_rate_term_sigma.csv', index_col=0).values
    results = pd.read_csv('european/npv_delta_gamma_vega_theta_rho.csv', index_col=0).values

    filtered_rows = ((np.log(params[:, 0]) < (params[:, 1] - 0.5 * params[:, 3]**2) * params[:, 2]
                      + 3 * params[:, 3] * np.sqrt(params[:, 2]))
                   & (np.log(params[:, 0]) > (params[:, 1] - 0.5 * params[:, 3]**2) * params[:, 2]
                      - 3 * params[:, 3] * np.sqrt(params[:, 2])))
    params = params[filtered_rows]
    results = results[filtered_rows]

    print(results.shape)

    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(params, results, test_size=0.1, random_state=2024)

    # 转换数据为张量
    X_train = torch.tensor(X_train.astype(np.float32)).to(device)
    y_train = torch.tensor(y_train.astype(np.float32)[:, [0]]).to(device)

    X_test = torch.tensor(X_test.astype(np.float32), requires_grad=True).to(device)
    y_test = torch.tensor(y_test.astype(np.float32)).to(device)

    # 构建TensorDataset和DataLoader
    train_dataset = TensorDataset(X_train, y_train)
    test_dataset = TensorDataset(X_test, y_test)

    batch_size = 64
    train_loader = DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True)
    test_loader = DataLoader(dataset=test_dataset, batch_size=1, shuffle=False)

    # 初始化模型、损失函数和优化器
    model = OptionPricingNN().to(device)
    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=0.001)

    # print("Training begins: ")
    #
    # # 训练模型
    # epochs = 10
    # for epoch in range(epochs):
    #     model.train()  # 设置模型为训练模式
    #     running_loss = 0.0
    #     for batch_idx, (data, target) in enumerate(train_loader):
    #         optimizer.zero_grad()
    #         output = model(data)
    #         loss = criterion(output, target)
    #         loss.backward()
    #         optimizer.step()
    #         running_loss += loss.item()
    #
    #     print(f'Epoch {epoch + 1}/{epochs}, Loss: {running_loss / len(train_loader)}')

    # torch.save(model.state_dict(), 'european/model_weights_4_50_50_1.pth')

    model.load_state_dict(torch.load('european/model_weights_4_50_50_1.pth', map_location=torch.device(device), weights_only=True))

    # 在测试集上评估模型
    print("Evaluation begins: ")
    start_time = time.time()  # 记录开始时间

    model.eval()
    criterion2 = DiffLoss()

    npv_loss, delta_loss, vega_loss, theta_loss, rho_loss = 0.0, 0.0, 0.0, 0.0, 0.0
    counter = 0

    for data, target in test_loader:
        if counter % 10000 == 0:
            print(counter)
        counter = counter + 1

        data, target = data.to(device), target.to(device)

        x1 = model.fc1(data).to(device)
        x1.retain_grad()
        output = model.layers(x1)

        loss = criterion2(output, target[:, [0]])
        npv_loss += loss.item()**2

        model.zero_grad()
        loss.backward()

        first_derivatives = torch.mm(x1.grad, model.fc1.weight)

        # 计算希腊字母
        delta = -data[0, 0] * first_derivatives[0, 0] + output[0, 0]
        vega = first_derivatives[0, 3]
        theta = -first_derivatives[0, 2]
        rho = first_derivatives[0, 1]

        target_delta = target[0, 1]
        target_vega = target[0, 3]
        target_theta = target[0, 4]
        target_rho = target[0, 5]

        delta_loss = delta_loss + (delta - target_delta) ** 2
        vega_loss = vega_loss + (vega - target_vega) ** 2
        theta_loss = theta_loss + (theta - target_theta) ** 2
        rho_loss = rho_loss + (rho - target_rho) ** 2

    print(f'NPV   Loss: {npv_loss / len(test_dataset)}')
    print(f'Delta Loss: {delta_loss / len(test_dataset)}')
    print(f'Vega  Loss: {vega_loss / len(test_dataset)}')
    print(f'Theta Loss: {theta_loss / len(test_dataset)}')
    print(f'Rho   Loss: {rho_loss / len(test_dataset)}')

    end_time = time.time()  # 记录结束时间
    print("执行时间：", end_time - start_time)

2. 局部波动率模型

实际当中，波动率是一个$T$和行权价$K$的函数，即所谓的隐含波动率曲面$\sigma (T,K)$，由此推出局部波动率曲面${\sigma }_{LV}(T,K)$，假设标的价格满足随机微分方程

$dS=rSdt+{\sigma }_{LV}(t,S)Sd{W}_t$

这就是局部波动率模型，这玩意儿的特点就是模拟到任何一个期限产生的价格分布与隐含波动率曲面所蕴含的风险中性下的概率分布(转化为欧式期权价格，用Butter Fly组合逼近)是一致的。

2.1 模型

对于欧式期权来说，我们可以直接将隐含波动率$\sigma (T,K)$代入Black-Scholes公式进行计算，在这里我们仍然将其归为局部波动率模型，并不是说我要给欧式期权上MC或PDE，最后跑出一个跟解析解一致的玩意儿沾沾自喜(但那是一个很好的练习)。对于路径依赖强的奇异期权，我们就不得不那么做了，但是输入还是$\sigma (T,K)$，${\sigma }_{LV}(t,S)$则是on-the-fly式的生成的，我假设你已经知道如何从隐含波动率曲面推出局部波动率曲面。

对于隐含波动率曲面的生成，经典的做法是将场内期权的报价转化成隐含波动率，对离散的点进行拟合或插值生成完整的波动率曲面。如果采用神经网络，我们恐怕得将上节中的$\sigma$扩充为一组刻画曲面的数字。我们分权益类和外汇类期权两种分别考虑。

外汇类期权市场数据，每一个点都是某个特定期限$T$和$\Delta$(这个遵循惯例的，跟你想要的那个可能不一样)上的隐含波动率值$\sigma (T,\Delta )$，任意时刻的报价期限和$\Delta$都是固定的，因此可以将二维的矩阵展成一维，与其他输入参数一起连成一个输入向量参与模型训练，对任何一个单点$\sigma (T,\Delta )$测算敏感度，我们就得到了Bucket Vega，天选之子。

权益类期权市场数据，每一个点都是某个特定期限$T$和行权价$K$上的看涨期权和看跌期权报价构成，在处理当中会结合远期价和虚值期权价格反推出$\sigma (T,K)$，但是，$(T,K)$的数目和分布并不规则。我们可以借鉴外汇期权的情形进行预处理，基于$\sigma (T,K)$先构造出完整的隐含波动率曲面(通常是先在$K$维度上进行样条插值，SVI或SABR拟合，然后在$T$维度上对${\sigma }^{2}T$进行线性插值)，然后获取一个固定网格上(由一组固定期限和一组固定行权价构成)的隐含波动率矩阵再进行训练。这里会有几个值得考虑的问题，

• 计算时(比如我们为了准备训练目标时做的)，在规范网格上二次插值产生的隐含波动率大概率无法复现市场数据值(选取一个较密的规范网格可以减小这种误差，但是会需要更大的计算量)，这是一个值得评估的模型误差。
• 如果我们在$K$维度上选择的是SVI或SABR拟合，那么每一条波动率微笑曲线的决定参数只有5个或更少，因此完全可以考虑将曲线拟合放在数据预处理环节，不同期限上的曲线参数放入模型输入向量。
• 前面的困扰其实都来自于我们对于训练数据的准备上，但是如果我们有了海量的历史数据，我们是不是连波动面插值都不需要了，期权参数也不需要自己预设了，直接拿来训练。外汇可能可以(当然，别在欧式期权上浪费时间了)，但权益类怎么办？我目前没想到。

我们下面遵循上面的第二个思路处理权益类期权。

由于输入向量的大幅延长，在目标值的准备上，需要的计算量都要大的多，但是这种大规模的计算储备对于日后的应用来说是划得来的。做个比方，这好比我们尽可能多的算出了各种市场参数搭配期权参数下的估值，我们把这些数据存了下来并进行了一种拟合(注意：这绝不是插值)，这种拟合可能要比我们直接进行高维空间插值更平滑和接近模型内在，最后，这种拟合要比常规的数值算法(MC或PDE)算起来更快。

2.2 SVI拟合

SVI(Stochastic Volatility Inspired)是由Gatheral提出的一种波动率微笑拟合方式，在权益期权上用的比较多。在有市场报价的期限$T$上，令$k=ln(K/F)$($F$表示远期价)，我们称之为$log-foward-moneyness$，拟合出波动率微笑曲线

${\sigma }^2(T,k)=a+b[\rho (k-m)+\sqrt{(k-m{)}^2+{\sigma }^2}]$

也有些人会把$T{\sigma }^2(T,k)$作为拟合对象，可以看出，它一共有5个参数$[a,b,\rho ,m,\sigma ]$，这5个参数可以转换成更有直觉意义的平值波动率，平值$Skew$，平值$Curvature$，左极限$Skew$，右极限$Skew$(有兴趣自己推导去)。当市场上有多个报价期限的时候，我们拟合出多组参数。为了避免$K$方向上的无风险套利，参数之间满足下列关系，

$b>0,$

$\sigma \ge 0,$

$|\rho |\le 1,$

$a\ge -b\sigma \sqrt{1-{\rho }^2},$

$b\le \frac{4}{(1+|\rho |)T}.$

我们现在假设市场上只有一个期限(这时候已经无所谓它是哪个期限上的了)，这样做的目的主要还是出于训练效率，并且，在点抽样的时候我们也不会受到期限无风险套利带来的参数限制。再简化，我们令$a=0.01,b=1.0,\sigma =0.05$，令无风险利率$r=0$。现在，我们的模型输入向量变成了$[K,T,\rho ,m]$。我们生成均匀的空间网格$linspace(0.5,1.5,101)\times linspace(0.01,1.01,101)\times linspace(-0.7,0.3,51)\times linspace(-0.3,0.3,61)$，这样，总样本点数大概3000万，用C++解析地算出$V$，${\Delta }_{Sticky-Moneyess}$，${\Delta }_{Sticky-Strike}$(定义见下文)，和$\Theta$，并存下来。随机选择其中的90%为训练集，余下的为测试集。

考虑一个欧式期权，当市场数据$S$发生变化，它的$moneyness$($K/S$)一定是变化的，$k$也是变化的，用变化的$k$从上面的式子中去索取隐含波动率，就是Sticky-Moneyness的观点：我只认$moneynees$。我们有时也习惯于将隐含波动率写成$\sigma (T,K)$，期权的$K$并没有变，那隐含波动率就不变，该吃吃，该喝喝，这就是Sticky-Strike的观点，这种观点下的$\Delta$和$\Gamma$就是我们简单的将隐含波动率代入Black-Scholes公式下得到的。

事实上，当我们采用Sticky-Strike的时候，SVI的参数似乎要重新拟合了，因为市场上报价的每个$k$都变了，但如果“货也要，钱又不想给”怎么办？好办，所有的$k$都发生了相同方向和大小的偏移，那么波动率微笑曲线无非就是发生一点左右平移，这可以通过调整SVI的参数$m$来实现。让我们来计算一下两种观点之下的$\Delta$，仍然假设$S$从$1$变化到$1+ϵ$，

${\Delta }_{Sticky-Moneyness}={\lim }_{ϵ\to 0}\frac{(1+ϵ)\hat{V}(K/(1+ϵ),m)-\hat{V}(K,m)}{ϵ}=-K\cdot \frac{\partial \hat{V}(K,m)}{\partial K}+\hat{V}(K,m),$

${\Delta }_{Sticky-Strike}={\lim }_{ϵ\to 0}\frac{(1+ϵ)\hat{V}(K/(1+ϵ),m-ln(1+ϵ))-\hat{V}(K,m)}{ϵ}=-K\cdot \frac{\partial \hat{V}(K,m)}{\partial K}+\hat{V}(K,m)-\frac{\partial \hat{V}(K,m)}{\partial m},$

$k$变化了$-ln(1+ϵ)$，为了取到不变的隐含波动率满足Sticky-Strike要求，$m$就得发生相应的调整。可见，两种$\Delta$都可以用AAD的方法来计算了。

2.3 数值

选择的神经网络大小为$4\times 50\times 50\times 1$，学习率设置为0.001，在训练集上前后传播10次。结果如下

	$V$	${\Delta }_{Sticky-Strike}$	${\Delta }_{Sticky-Moneyness}$	$\Theta$
$MSE$	$6.25\times 10^{-7}$	$5.51\times 10^{-4}$	$8.92\times 10^{-4}$	$3.75\times 10^{-4}$

算的稍微有点慢(我用的Air)，连训练带测试大约耗时24小时。欣慰的是，结果算是挺完美了。

2.4 代码

老规矩，csv文件自己用C++去准备。

european_call_svi.py

import math
import time
import scipy.stats as stats

import numpy as np
from sklearn.model_selection import train_test_split
from torch import nn, optim
import torch
from torch.utils.data import TensorDataset, DataLoader
import pandas as pd

# 设置设备
device = torch.device('mps' if torch.backends.mps.is_available() else 'cpu')


def black_scholes_call(K, T, rho, m):
    S = 1
    r = 0
    a = 0.01
    b = 1.0
    sigma = 0.05

    F = S * math.exp(r*T)
    k = math.log(K / F)
    x = k - m

    var = a + b * (rho * x + math.sqrt(x * x + sigma * sigma))
    vol = math.sqrt(var)

    d1 = (math.log(S / K) + (r + 0.5 * var) * T) / (vol * math.sqrt(T))
    d2 = d1 - vol * math.sqrt(T)

    call_price = S * stats.norm.cdf(d1) - K * math.exp(-r * T) * stats.norm.cdf(d2)
    delta_ss = stats.norm.cdf(d1)
    vega = S * stats.norm.pdf(d1) * math.sqrt(T) / math.sqrt(2 * math.pi)
    theta = -S * vol * vol / (2 * math.sqrt(T)) - r * K * math.exp(-r * T) * stats.norm.cdf(d2)

    dvar_dk = b * (rho + x / math.sqrt(x * x + sigma * sigma))
    dvol_dk = dvar_dk / (2.0 * vol)
    delta_sm = delta_ss - vega * dvol_dk / S

    return call_price, delta_ss, delta_sm, theta


# 构建神经网络模型
class OptionPricingNN(nn.Module):
    def __init__(self, width=50):
        super(OptionPricingNN, self).__init__()
        self.width = width
        self.fc1 = nn.Linear(4, width)
        self.layers = nn.Sequential(
            nn.ReLU(),
            nn.Linear(width, width),
            nn.ReLU(),
            nn.Linear(width, 1)
        )

    def forward(self, x):
        x = self.fc1(x)
        x = self.layers(x)

        return x


class DiffLoss(nn.Module):
    def __init__(self):
        super(DiffLoss, self).__init__()

    def forward(self, input, target):
        # 计算MSE损失
        loss = input - target

        return loss


if __name__ == "__main__":
    params = pd.read_csv('european/strike_term_rho_m_svi.csv', index_col=0).values
    results = pd.read_csv('european/npv_ss_sm_theta_svi.csv', index_col=0).values

    print(results.shape)

    # 划分训练集和测试集
    X_train, X_test, y_train, y_test = train_test_split(params, results, test_size=0.1, random_state=2024)

    # 转换数据为张量
    X_train = torch.tensor(X_train.astype(np.float32)).to(device)
    y_train = torch.tensor(y_train.astype(np.float32)[:, [0]]).to(device)

    X_test = torch.tensor(X_test.astype(np.float32), requires_grad=True).to(device)
    y_test = torch.tensor(y_test.astype(np.float32)).to(device)

    # 构建TensorDataset和DataLoader
    train_dataset = TensorDataset(X_train, y_train)
    test_dataset = TensorDataset(X_test, y_test)

    batch_size = 64
    train_loader = DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True)
    test_loader = DataLoader(dataset=test_dataset, batch_size=1, shuffle=False)

    # 初始化模型、损失函数和优化器
    model = OptionPricingNN().to(device)
    criterion = nn.MSELoss()
    optimizer = optim.Adam(model.parameters(), lr=0.001)

    print("Training begins: ")

    epochs = 10
    for epoch in range(epochs):
        model.train()  # 设置模型为训练模式
        running_loss = 0.0
        for batch_idx, (data, target) in enumerate(train_loader):
            optimizer.zero_grad()
            output = model(data)
            loss = criterion(output, target)
            loss.backward()
            optimizer.step()
            running_loss += loss.item()

        print(f'Epoch {epoch + 1}/{epochs}, Loss: {running_loss / len(train_loader)}')

    torch.save(model.state_dict(), 'european/model_weights_4_50_50_1_svi.pth')

    # model.load_state_dict(torch.load('european/model_weights_4_50_50_1_svi.pth', map_location=torch.device(device), weights_only=True))

    # 在测试集上评估模型
    print("Evaluation begins: ")
    start_time = time.time()  # 记录开始时间

    model.eval()
    criterion2 = DiffLoss()

    npv_loss, delta_ss_loss, delta_sm_loss, theta_loss = 0.0, 0.0, 0.0, 0.0
    counter = 0

    for data, target in test_loader:
        if counter % 10000 == 0:
            print(counter)
        counter = counter + 1

        data, target = data.to(device), target.to(device)

        x1 = model.fc1(data).to(device)
        x1.retain_grad()
        output = model.layers(x1)

        loss = criterion2(output, target[:, [0]])
        npv_loss += loss.item()**2

        model.zero_grad()
        loss.backward()

        first_derivatives = torch.mm(x1.grad, model.fc1.weight)

        # 计算希腊字母
        delta_sm = -data[0, 0] * first_derivatives[0, 0] + output[0, 0]
        delta_ss = delta_sm - first_derivatives[0, 3]
        theta = -first_derivatives[0, 1]

        target_delta_ss = target[0, 1]
        target_delta_sm = target[0, 2]
        target_theta = target[0, 3]

        delta_ss_loss = delta_ss_loss + (delta_ss - target_delta_ss) ** 2
        delta_sm_loss = delta_sm_loss + (delta_sm - target_delta_sm) ** 2
        theta_loss = theta_loss + (theta - target_theta) ** 2

    print(f'NPV     Loss: {npv_loss / len(test_dataset)}')
    print(f'DeltaSS Loss: {delta_ss_loss / len(test_dataset)}')
    print(f'DeltaSM Loss: {delta_sm_loss / len(test_dataset)}')
    print(f'Theta   Loss: {theta_loss / len(test_dataset)}')

    end_time = time.time()  # 记录结束时间
    print("执行时间：", end_time - start_time)

3. 总结

我们将简单神经网络应用在了欧式期权的估值上，分别基于Black-Scholes模型和局部波动率模型，采用AAD的思想来计算一些希腊字母，这是个很有意思的尝试，也让我们似乎看到了一线在复杂结构上应用的曙光。

慢着，还是得想清楚。首先，我们该如何规范市场行情和合约参数来作为输入。关于市场行情，前面提到了权益类的隐含波动率曲面就是个麻烦。至于合约参数，那会更让人挠头，了解一点具体的场外期权合约就知道了。其次，训练数据的准备，市场报价噪声太多，我们恐怕还是得调基础库来算，这个计算量不会小，有些产品没什么交易量，不值得去为之训练。用GPU来跑蒙特卡洛可能是一个不错的选择，可以先让ChatGPT写一个简单的例子再慢慢攒。总之，别为了神经网络发神经。要是前面的问题都解决了，那就像调鸡尾酒一样调调模型参数吧。我目前能想到的最适合深度学习来训练的对象应该是外汇的障碍和触碰期权，那玩意儿结构相对固定，交易量大，讲究一点的都会用局部随机波动率(Local Stochastic Volatility Model)模型来做，而那是一个恶心的双因子模型，计算相当耗时。

关于ChatGPT，确实是一个很有用的工具，有一些想法可以立马去尝试，以前我们很多人是面向百度编程，现在得益于技术的发展升级了。但最后的结果还是靠人为干涉，这世界怎么能没有可爱的程序员呢？