Black-Scholes Equation in Python 改良版

前回のコードをある程度改良しましたので、載せておきたいと思います。主な改良点は以下。

• ヨーロピアン(コール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)

て、そんなことやってる場合じゃない。早くテスト対策行なわないと。