Golang 期权 指标计算

用golang做期权行情指标计算时，关于sigma, delta, gamma, vega, theta, rho等的计算，需要用到black-scholes模型

期权指标计算其他语言已有现成的包或者代码，我参考了其他博主的python代码，用go完整的写了一遍，包含隐含波动率、正态分布函数等

 直接贴代码

//DN 标准正态分布的概率密度函数公式
func DN(x float64) float64 {
	return math.Exp(-0.5*math.Pow(x, 2)) / math.Sqrt(2*math.Pi)
}

//CalcVega 计算期权Vega指标  隐含波动率变化引起的期权价格变化
/*
参考:  https://blog.youkuaiyun.com/xiaowu1997/article/details/122279216
St 	  标的资产在t时刻的价格-标的最新价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
*/
func CalcVega(St, K, T, r, sigma float64) float64 {
	d1 := calcD1(St, K, sigma, r, T)
	return calcVega(d1, St, T)
}

func calcVega(d1, St, T float64) float64 {
	return St * math.Sqrt(T) * DN(d1) / 100.0
}

//CalcDelta 标的资产价格变化引起期权价格变化 计算期权Delta值  n=1看涨期权   n=-1看跌期权
/*
St 	  标的资产在t时刻的价格-标的最新价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
n     n=1 看涨  n=-1看跌
*/
func CalcDelta(St, K, T, r, sigma float64, n float64) float64 {
	d1 := calcD1(St, K, sigma, r, T)
	return calcDelta(d1, n)
}

func calcDelta(d1, n float64) float64 {
	return n * N(n*d1)
}

//CalcGamma 标的资产价格变化引起delta值的变化
/*
St 	  标的资产在t时刻的价格-标的最新价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
*/
func CalcGamma(St, K, T, r, sigma float64) float64 {
	d1 := calcD1(St, K, sigma, r, T)
	return calcGamma(d1, St, T, sigma)
}

/*
标的资产价格变化引起delta值的变化
d1 		期权对股价的敏感程度
St 	  	标的资产在t时刻的价格-标的最新价
T	  	有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
sigma 	隐含波动率
*/
func calcGamma(d1, St, T, sigma float64) float64 {
	return DN(d1) * (St * sigma * math.Sqrt(T)) / 1000
}

//CalcRho 期权价格对（无风险）利率变化的敏感程度
/*
参考:  https://blog.youkuaiyun.com/xiaowu1997/article/details/122347991
St 	  标的资产在t时刻的价格-标的最新价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
n     n=1 看涨  n=-1看跌
*/
func CalcRho(St, K, T, r, sigma, n float64) float64 {
	d1 := calcD1(St, K, sigma, r, T)
	d2 := calcD2(d1, sigma, T)
	return calcRho(d2, K, T, r, n)
}

/*
标的资产价格变化引起delta值的变化
d2	期权最后被执行的可能性
K	执行价
T	有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r	无风险连续复利利率
n	n=1 看涨  n=-1看跌
*/
func calcRho(d2, K, T, r, n float64) float64 {
	return n * K * T * math.Exp(-r*T) * N(n*d2) / 100
}

//CalcTheta 期权的时间价值随时间流逝耗损的速度
/*
参考 https://nbviewer.org/github/poppinkingone/jupyter_test/blob/master/Option.ipynb#put-:delta-=-N(-d1)
St 	  标的资产在t时刻的价格-标的最新价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
n     n=1 看涨  n=-1看跌
*/
func CalcTheta(St, K, T, r, sigma float64, n float64) float64 {
	d1 := calcD1(St, K, sigma, r, T)
	d2 := calcD2(d1, sigma, T)
	//业内习惯Theta用来衡量每日的time decay，而BS Model中的时间单位是年，所以按此公式算出来的Theta需要再除以365
	return calcTheta(d1, d2, St, K, T, r, sigma, n)
}

/*
期权的时间价值随时间流逝耗损的速度
d1 		期权对股价的敏感程度
d2		期权最后被执行的可能性
St 	  	标的资产在t时刻的价格-标的最新价
K		执行价
T		有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r		无风险连续复利利率
sigma 	隐含波动率
n		n=1 看涨  n=-1看跌
*/
func calcTheta(d1, d2, St, K, T, r, sigma float64, n float64) float64 {
	//业内习惯Theta用来衡量每日的time decay，而BS Model中的时间单位是年，所以按此公式算出来的Theta需要再除以365
	return (-St*DN(d1)*sigma/(2*math.Sqrt(T)) - n*r*K*math.Exp(-r*T)*N(n*d2)) / 365
}

// N 计算值u的标准正态分布概率
func N(u float64) float64 {
	//取绝对值
	y := math.Abs(u)
	//平方值
	y2 := math.Pow(y, 2)
	z := math.Exp(-0.5*y2) * 0.398942280401432678
	p := float64(0)
	k := 28
	s := float64(-1)
	fj := float64(k)
	if y > 3 {
		//当y>3
		for i := 1; i <= k; i++ {
			p = fj / (y + p)
			fj = fj - 1.0
		}
		p = z / (y + p)
	} else {
		//当y<3
		for i := 1; i <= k; i++ {
			p = fj * y2 / (2.0*fj + 1.0 + s*p)
			s = -s
			fj = fj - 1.0
		}
		p = 0.5 - z*y/(1-p)
	}
	if u > 0 {
		p = 1.0 - p
	}
	return p
}

//计算期权对股价的敏感程度 BS mode中的d1
/*
S     标的资产在T时刻的价格
K	  执行价
sigma 隐含波动率
r 	  无风险连续复利利率
t	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
*/
func calcD1(S float64, K float64, sigma float64, r float64, t float64) float64 {
	return (math.Log(S/K) + (r+0.5*math.Pow(sigma, 2))*t) / (sigma * math.Sqrt(t))
}

//计算期权最后被执行的可能性 BS mode中的d2
/*
d1 		期权对股价的敏感程度
sigma 	隐含波动率
t  		期权有效天数与一年365天的比值
*/
func calcD2(d1 float64, sigma float64, t float64) float64 {
	return d1 - sigma*math.Sqrt(t)
}

// GetC 计算预期定价
/*
S 		标的资产在T时刻的价格
K 		期权交割价格
sigma	一年度化方差
r 		连续复利无风险利率
t 		期权有效期
*/
func GetC(S float64, K float64, sigma float64, r float64, t float64) float64 {
	d1 := calcD1(S, K, sigma, r, t)
	d2 := calcD2(d1, sigma, t)
	return S*N(d1) - K*math.Exp(-r*t)*N(d2)
}

// CalcVolatility 计算引伸波幅,对冲值。采取逼近法计算
/*
isRise 是看涨还是看跌 认购true； 认沽false）
warPx 期权现价
cr    换股比率
S     标的资产在T时刻的价格
r 	  无风险连续复利利率
L	  实施价格(行权价格)
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
Dt    期权的期内预期分红（注意分红有特定时间的，等不到分红的期权预期分红就是零）
*/
func CalcVolatility(isRise bool, warPx float64, cr float64, S float64, r float64, L float64, T float64, Dt float64) (float64, float64) {
	if warPx == 0 || S == 0.0 {
		return 0.0, 0.0
	}
	NowC := warPx * cr
	sigma := float64(0)
	var c float64
	for i := 0; i < 20; i++ {
		c, _ = CalcC(isRise, S, r, sigma*sigma, L, T, Dt)
		if c > NowC {
			sigma -= 0.1
			break
		} else {
			sigma += 0.1
		}
	}
	for i := 0; i < 10; i++ {
		c, _ = CalcC(isRise, S, r, sigma*sigma, L, T, Dt)
		if c > NowC {
			sigma -= 0.01
			break
		} else {
			sigma += 0.01
		}
	}
	for i := 0; i < 10; i++ {
		c, _ = CalcC(isRise, S, r, sigma*sigma, L, T, Dt)
		if c > NowC {
			sigma -= 0.001
			break
		} else {
			sigma += 0.001
		}
	}
	var delta float64
	for i := 0; i < 10; i++ {
		c, delta = CalcC(isRise, S, r, sigma*sigma, L, T, Dt)
		if c > NowC {
			break
		} else {
			sigma += 0.0001
		}
	}
	return sigma, delta
}

//CalcC 计算理论价
/*
isRise true看涨  false 看跌
S      标的资产在T时刻的价格
r 	   无风险连续复利利率
sigma  隐含波动率
K	   执行价
T	   有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
Dt     期权的期内预期分红（注意分红有特定时间的，等不到分红的期权预期分红就是零）
*/
func CalcC(isRise bool, S float64, r float64, sigma float64, K float64, T float64, Dt float64) (float64, float64) {
	d1 := CalcD1(S, r, sigma, K, T)
	d2 := CalcD2(d1, sigma, T)
	a := Dt
	b := N(d1)
	c := math.Exp(-1 * r * T)
	d := N(d2)
	C := (S-a)*b - K*c*d
	if !isRise {
		C = K*c*(1-d) - (S-a)*(1-b)
		b = b - 1
	}
	return C, b
}

//CalcD1 计算期权对股价的敏感程度 BS mode中的d1
/*
S     标的资产在T时刻的价格
r 	  无风险连续复利利率
sigma 隐含波动率
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
*/
func CalcD1(S float64, r float64, sigma float64, K float64, T float64) float64 {
	pv := K / math.Pow(1+r, T)
	return math.Log(S/pv)/(math.Sqrt(sigma*T)) + math.Sqrt(sigma*T)/2
}

//CalcD2 计算期权最后被执行的可能性 BS mode中的d2
/*
d1 		期权对股价的敏感程度
sigma 	隐含波动率
T  		期权有效天数与一年365天的比值
*/
func CalcD2(d1 float64, sigma float64, T float64) float64 {
	return d1 - math.Sqrt(sigma*T)
}

/*
计算看涨期权理论价格
S     标的资产在T时刻的价格
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
*/
func bsCall(S, K, T, r, sigma float64) float64 {
	d1 := (math.Log(S/K) + (r+0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	d2 := (math.Log(S/K) + (r-0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	return S*N(d1) - K*math.Exp(-r*T)*N(d2)
}

/*
计算期权理论价格
S     标的资产在T时刻的价格
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
n     n=1 看涨  n=-1看跌
*/
func bsPrice(S, K, T, r, sigma, n float64) float64 {
	d1 := (math.Log(S/K) + (r+0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	d2 := (math.Log(S/K) + (r-0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	return n*S*N(n*d1) - n*K*math.Exp(-r*T)*N(n*d2)
}

/*
计算看跌期权理论价格
S     股票现价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
sigma 隐含波动率
*/
func bsPut(S, K, T, r, sigma float64) float64 {
	d1 := (math.Log(S/K) + (r+0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	d2 := (math.Log(S/K) + (r-0.5*math.Pow(sigma, 2))*T) / (sigma * math.Sqrt(T))
	return -S*N(-d1) + K*math.Exp(-r*T)*N(-d2)
}

//CalcSigma 计算隐含波动率
/*
参考1: https://blog.youkuaiyun.com/hzk427/article/details/104501847/
参考2: https://blog.youkuaiyun.com/qq_18822147/article/details/108037407?utm_medium=distribute.pc_relevant.none-task-blog-2~default~baidujs_baidulandingword~default-0-108037407-blog-104501847.pc_relevant_multi_platform_whitelistv1_exp2&spm=1001.2101.3001.4242.1&utm_relevant_index=3
S     股票现价
K	  执行价
T	  有效期，期权有效天数与一年365天的比值。如果期权有效期为100天，则T=100/365=0.274
r 	  无风险连续复利利率
price 期权现价
n     n=1 看涨，  n=-1看跌
*/
func CalcSigma(S, K, T, r, price, n float64) float64 {
	EPS := 1e-5
	lo := 1e-5
	hi := 0.5
	var c float64
	for lo+EPS < hi {
		sigma := (lo + hi) / 2
		if n == 1 {
			c = bsCall(S, K, T, r, sigma)
		} else {
			c = bsPut(S, K, T, r, sigma)
		}
		//c = bsPrice(S, K, T, r, sigma, n)
		if c-price < 0 {
			lo = sigma
		} else {
			hi = sigma
		}
	}
	return lo
}

//CalcOptionQuota 计算期权指标
/*
参考：https://blog.youkuaiyun.com/Zita_11/article/details/104200887
S 		标的股票现价
K 		期权执行价
T 		剩余到期日比例
r 		无风险连续复利利率
price 	期权现价
isRise 	看涨true/看跌false
*/
func CalcOptionQuota(S, K, T, r, price float64, isRise bool) (sigma, delta, gamma, vega, theta, rho float64) {
	n := 1.0
	if !isRise {
		n = -1.0
	}
	sigma = CalcSigma(S, K, T, r, price, n)
	d1 := calcD1(S, K, sigma, r, T)
	d2 := calcD2(d1, sigma, T)
	delta = calcDelta(d1, n)
	gamma = calcGamma(d1, S, T, sigma)
	vega = calcVega(d1, S, T)
	theta = (-S*DN(d1)*sigma/(2*math.Sqrt(T)) - n*r*K*math.Exp(-r*T)*N(n*d2)) / 365
	rho = calcRho(d2, K, T, r, n)
	return
}