python / scipy 中的牛顿法

Question

我正在尝试在Python上使用牛顿算法来获取函数的根。即使我更改精度级别，我也会出现运行时错误。您能帮我了解如何改进吗？

最好的， 国标

在我的“简单”代码和根查找部分下面：

from scipy.stats import norm
import numpy as np 
import matplotlib.pyplot as plt  
import scipy.optimize as opt


def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

x0=112
r=0
T=1
d=0
sigma=0.2
q=0
K1=110
K2=130




res2= opt.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)


the error i get is: res2= opt.zeros.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)
/Users/anaconda/lib/python3.5/site-packages/scipy/optimize/zeros.py:173: RuntimeWarning: Tolerance of 0.011300000000005639 reached
  warnings.warn(msg, RuntimeWarning)

Answer 1

在如此稀疏的背景下很难提供建议。但一些评论：

答：

您没有使用牛顿，如此处所述：

如果提供了func的导数fprime，则使用Newton-Raphson方法，否则使用割线方法。

乙：

你的错误来自这里我会说：

由于这些 tol 值是硬编码的，我想这不应该发生！

# Secant method
p0 = x0
if x0 >= 0:
    p1 = x0*(1 + 1e-4) + 1e-4
else:
    p1 = x0*(1 + 1e-4) - 1e-4
q0 = func(*((p0,) + args))
q1 = func(*((p1,) + args))
for iter in range(maxiter):
    if q1 == q0:
        if p1 != p0:
            msg = "Tolerance of %s reached" % (p1 - p0)
            warnings.warn(msg, RuntimeWarning)
        return (p1 + p0)/2.0

含义：您的代码可能有问题！

C：

让我们尝试更慢但更安全的二分法：

# brackets not tuned! it's just some trial!
res2 = opt.bisect(vega_callspread, 0, 200, args=(r,T,d,sigma,q,K1,K2))

输出：

X:\so_newton.py:9: RuntimeWarning: divide by zero encountered in log
  d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
X:\so_newton.py:10: RuntimeWarning: divide by zero encountered in log
  d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

这是一个坏兆头。

D（跟随C）：

您的功能如下：

def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

然后你打电话：

x0=112
r=0
T=1
d=0
sigma=0.2
q=0
K1=110
K2=130

args=(r,T,d,sigma,q,K1,K2,)

没有S！

因此，要么您危险地重命名变量，要么您这边出现了一些错误！

S

的值是

log(S/K1)

和co中的问题。

这不是关于

S/K1

而是关于这个：

import numpy as np
np.log(0)
# __main__:1: RuntimeWarning: divide by zero encountered in log
# -inf

我从这个SO答案得到的。

您现在可以尝试解释这对您的函数做了什么，因为此 log-eval 将获得

-np.inf

的值。

E scipy 如何处理这个问题（遵循 D）：

我懒得阅读 args-handling 的文档/源代码，但让我们检查一些东西（向 func 添加打印；像之前一样使用二分法）：

def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    print('S: ', S)

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

输出：

S:  0
X:\so_newton.py:11: RuntimeWarning: divide by zero encountered in log
  d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
X:\so_newton.py:12: RuntimeWarning: divide by zero encountered in log
  d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
S:  0

所以看来 arg 处理是通过 arg-name 进行的（而不仅仅是调用中的参数位置；我在这里缺少正确的编程语言术语；又懒了！）

这（S=0）可能非常糟糕（不是你想要的），这是你的错误！

编辑

发表评论后，看来您正在尝试优化

S

。这让我很清楚，您正在使用一些优化算法来优化

x

，而您的函数中没有

x

！

我不是在这里分析你的任务，但你可能想让你的函数使用一些

x

（由

x0

初始化），因为这是

scipy.optimize

的总体思路。我们可以保留名称

S

，但它需要是函数的第一个参数。这一切都在文档中进行了解释。

所以：

def vega_callspread(S, r,T,d,sigma,q,K1,K2):  # S now first argument !!!
    ...

res2= opt.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)  # S removed from args; S is init by x0 -> read docs!

输出：

117.214682594

Answer 2

S

是常数，将其作为常数（每个特定分析日期的库存 Pr），除以

K

（罢工 Pr） - 不，您不会乘以零（如标记的答案中所示）， - 来自该领域的本质（期权定价）。尝试更改您的容差参数，例如

tol=1e-8