前回のコードをある程度改良しましたので、載せておきたいと思います。主な改良点は以下。
- ヨーロピアン(コールorプット)オプションのデルタ、ガンマ、ベガ、シータ、ローを求める関数を追加しました。
- ヨーロピアン(コールorプット)オプションのインプライド・ボラティリティをニュートン法より求める関数を追加しました。
- その他若干の名称変更や、引数に関する説明を加えました。
数式は"Black-Scholes (1973) Option Pricing Formula"を参考にしました。
インプライド・ボラティリティに関しては適切な引数を加えないと、負の値や無限大に発散することがあるので注意が必要です。とはいえ実際の市場価格と理論価格がそこまで乖離することもないでしょうから、大抵正常な値を返すはずです。
商用、非商用に関わらず自由に使用、改変していいですし連絡する必要もないですが、無保証です。使う人はあんまいないと思いますが、一応建前として。
#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
Functions for calculating the Black-Scholes equation.
"""
import math
A1 = 0.319381530
A2 = -0.35653782
A3 = 1.781477937
A4 = -1.821255978
A5 = 1.330274429
RSQRT2PI = 0.3989422804
GAMMA = 0.2316419
def BS_Call(S,K,T,R,V):
"""Calculate the Black-Scholes formula for a call option.
Arguments:
S:the price of the underlying stock
K:the strike price
R:the continuously compounded risk free interest rate
T:the time in years until the expiration of the option
V:the volatility for the underlying stock
"""
d1,d2 = _getD1D2(S,K,T,R,V)
ret = S*_normalDist(d1)-K*math.exp(-R*T)*_normalDist(d2)
return ret
def BS_Put(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
ret = -S*_normalDist(-d1)+K*math.exp(-R*T)*_normalDist(-d2)
return ret
def BS_Call_IV(S,K,T,R,Pr,HV,maxNum=100,i=1):
"""Calculate the Implied Volatility using Newton's method.
Arguments:
Pr:the real price of call option
HV:the Historical Volatility
maxNum:the number of repeating calculation
i:index(Do not specify this parameter)
"""
if i > maxNum:
return HV
else:
newV = HV - (BS_Call(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
return BS_Call_IV(S,K,T,R,Pr,newV,maxNum,i+1)
def BS_Put_IV(S,K,T,R,Pr,HV,maxNum=10,i=1):
if i > maxNum:
return HV
else:
newV = HV - (BS_Put(S,K,T,R,HV)-Pr)/BS_Vega(S,K,T,R,HV)
return BS_Put_IV(S,K,T,R,Pr,newV,maxNum,i+1)
def BS_Call_Delta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1)
def BS_Put_Delta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1) - 1.0
def BS_Gamma(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return _normalDist(d1)/(S*V*math.sqrt(T))
def BS_Vega(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return S*_normalDist(d1)*math.sqrt(T)
def BS_Call_Theta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -S*_normalDist(d1)*V/(2*math.sqrt(T))-R*K*math.exp(-R*T)*_normalDist(d2)
def BS_Put_Theta(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -S*_normalDist(d1)*V/(2*math.sqrt(T))+R*K*math.exp(-R*T)*_normalDist(-d2)
def BS_Call_Rho(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return K*T*math.exp(-R*T)*_normalDist(d2)
def BS_Put_Rho(S,K,T,R,V):
d1,d2 = _getD1D2(S,K,T,R,V)
return -K*T*math.exp(-R*T)*_normalDist(-d2)
def _getD1D2(S,K,T,R,V):
vSqrtT = V*math.sqrt(T)
d1 = (math.log(float(S)/K)+(R+(0.5*V*V))*T)/(vSqrtT)
d2 = d1 - vSqrtT
return (d1,d2)
def _normalDist(x):
if x >= 0:
k = 1.0/(1.0 + GAMMA*x)
cnd = RSQRT2PI*math.exp(-0.5*x*x)*(k *(A1+ k*(A2+ k*(A3 + k*(A4 + k*A5)))))
ret = 1.0 - cnd
return ret
else:
return 1.0 - _normalDist(-x)
て、そんなことやってる場合じゃない。早くテスト対策行なわないと。
