局部波动率模型的快速蒙特卡洛

如何解决局部波动率模型的快速蒙特卡洛

我想通过蒙特卡洛模拟来计算期权价格。遵循SDE，实现的模型是马尔可夫过程：

d X_t = alpha * dt + beta^(1/2) * d W_t

其中alpha住在R ^ 2中，beta住在R ^（2 * 2）中，并且为正半定数。（^（1/2）-> cholesky）因此，X_t在R ^ 2中。

在使用网格网格上的有限差分（生成数组sigmasqX，sigmasqY和sigmaXY进行HJB方程解析之后，我获得了beta和alpha的值，然后将其插值为在(t,X_t)处获取beta和alpha值）。

我的问题是，我随后运行了蒙特卡洛仿真来获得期权价格，以得出隐含的波动率，但是这种仿真永远只能用于非常有限的20000路径仿真。

以下是我策略的主要内容：

我从首字母X_0
我生成一个大小为t_size*sample的普通样本，其中tsize是我的时间向量的长度，sample是生成的路径数
我在spicy.RectBivariateSpline数组上使用sigma进行插值。这是嵌套的for循环：

for j in np.arange(sample):

        X = np.array([np.log(var.S0),var.x2_0])
        broke = False

        for i in np.arange(1,ind[0]+2) :
            sigsqX  = L0[i-1]
            sigsqY  = L1[i-1]
            sigXY   = L2[i-1]
            
            beta[0,0] = sigsqX(X[0],X[1])
            beta[1,1] = sigsqY(X[0],X[1])
            beta[0,1] = sigXY(X[0],0] = beta[0,1]
            
            if not isPD(beta) :
                beta = nearestPD(beta)
                
            alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)

            bet_sq = np.linalg.cholesky(beta)
            X[0] =  X[0] + alpha[0]*dt + np.dot(bet_sq,W[j,i-1])[0]*np.sqrt(dt)
            X[1] =  X[1] + alpha[1]*dt + np.dot(bet_sq,i-1])[1]*np.sqrt(dt)
            
            if X[0] > x1vec[-1] or X[0] < x1vec[0]:
                broke = True
                break
                
            if X[1] > x2vec[-1] or X[1] < x2vec[0]:
                broke = True
                break

        if broke == False :
            vix_v_t1.append(X[1])
            
            for i in np.arange(ind[0]+2,ind[1]+2) :

                sigsqX  = L0[i-1]
                sigsqY  = L1[i-1]
                sigXY   = L2[i-1]
                
                beta[0,X[1])
                beta[1,X[1])
                beta[0,1]
                
                if not isPD(beta) :
                    beta = nearestPD(beta)
                    
                alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)
                bet_sq = np.linalg.cholesky(beta)
                X[0] =  X[0] + alpha[0]*dt + np.dot(bet_sq,i-1])[0]*np.sqrt(dt)
                X[1] =  X[1] + alpha[1]*dt + np.dot(bet_sq,i-1])[1]*np.sqrt(dt)

                if X[0] > x1vec[-1] or X[0] < x1vec[0]:
                    broke = True
                    break
                if X[1] > x2vec[-1] or X[1] < x2vec[0]:
                    broke = True
                    break

            if broke == False :
                spx_t1.append(X[0])
                for i in np.arange(ind[1]+2,ind[2]+2) :

                    sigsqX  = L0[i-1]
                    sigsqY  = L1[i-1]
                    sigXY   = L2[i-1]
                    
                    beta[0,X[1])
                    beta[1,X[1])
                    beta[0,1]
                    
                    if not isPD(beta) :
                        beta = nearestPD(beta)
                        
                    alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)
                    bet_sq = np.linalg.cholesky(beta)
                    X[0] =  X[0] + alpha[0]*dt + np.dot(bet_sq,i-1])[0]*np.sqrt(dt)
                    X[1] = X[1] + alpha[1]*dt + np.dot(bet_sq,i-1])[1]*np.sqrt(dt)
                    if X[0] > x1vec[-1] or X[0]<x1vec[0]:
                        broke = True
                        break
                    if X[1] > x2vec[-1] or X[1]<x2vec[0]:
                        broke = True
                        break

                if broke == False :
                    spx_t2.append(X[0])

    spot = [spx_t1,spx_t2,vix_v_t1]
    return spot

一些注意事项：

我检查是否beta，并根据需要对其进行调整，由于数值逼近，它有时具有很小的负特征值，并破坏了cholesky分解
有3个线性for循环可遍历整个时间集，因为选项具有三个到期日，因此更易于管理（此处无影响）
有时X_t的值会在网格之外徘徊，并破坏了模拟（由于外推误差），所以我实施了所有这些验证（if broke : break）
L0，L1和L2变量是python list类型的Python，它们在时间向量的每个日期保存网格的插值。

我知道为什么这段代码慢：

许多if语句使速度大大降低
为Beta查找最接近的PD可能会很昂贵
使用局部波动性方法，我无法根据NumPy的规则，蒙特卡洛过程进行并行化

我检查了一下，得出的插值微不足道（平均0.5秒）。但是，此路径仿真需要花费75秒才能完成2万条路径...

我的问题：有没有办法加快整个过程的速度？我无法想象我陷入了如此糟糕的执行时间...（与解决向后定价的PDE相比）以扩展计算能力而臭名昭著）

谢谢！

局部波动率模型的快速蒙特卡洛

如何解决局部波动率模型的快速蒙特卡洛

相关推荐