问题描述
我想通过蒙特卡洛模拟来计算期权价格。遵循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变量是pythonlist类型的Python，它们在时间向量的每个日期保存网格的插值。
我知道为什么这段代码慢：
- 许多
if语句使速度大大降低 - 为Beta查找最接近的PD可能会很昂贵
- 使用局部波动性方法，我无法根据NumPy的规则，蒙特卡洛过程进行并行化
我检查了一下，得出的插值微不足道（平均0.5秒）。但是，此路径仿真需要花费75秒才能完成2万条路径...
我的问题：有没有办法加快整个过程的速度？我无法想象我陷入了如此糟糕的执行时间...（与解决向后定价的PDE相比）以扩展计算能力而臭名昭著）
谢谢！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）