使用牛顿法和二分法计算期权的隐含波动率，把两种算法写在dolphindb的modules的opt_implied_volatility模块里。

2022.06.28日收盘,计算50etf执行价为3.0的看涨期权隐含波动率18.3。牛顿法耗时0.589ms,二分法耗时4.743ms

use opt_implied_volatility
optprice = 0.085//执行价为3，50etf看涨期权价格
K = 3.//行权价
T = (2022.07.27-2022.06.29) \365.
S = 3.043//2022.06.28日50etf收盘价格
r = 0.0002//无风险利率
opt_type = 'C'  // 看涨期权
timer implied_vol = implied_volatility_Newton(optprice, opt_type, S, K, T, r)/////Time elapsed: 0.589 ms
print(round(implied_vol * 100,2))

timer implied_vol=implied_volatility_Binary(optprice, S, K, r, T, opt_type)////Time elapsed: 4.743 ms
print(round(implied_vol * 100,2))

print(optprice)
print(round(bs_price(opt_type, S, K, T, r, implied_vol),3))

opt_implied_volatility模块如下

module opt_implied_volatility
def pdfNorm(x,mean=0,var=1){
	return exp(-square(x -mean) / (2 * square(var))) / (sqrt(2 * pi) * var)
 }
def bs_price(options_type, spot, strike, maturity, r, vol, q=0.0){
	/*
	 * options_type:'期权类型 ，C' or 'P'
	 * spot：现货价，float:
	 * strike：行权价，float 
	 * maturity：剩余天数，float、int
	 * r：无风险利率，float
	 * vol：波动率，float
	*/
	 if (options_type == 'C'){
	 	phi = 1}
	 else if(options_type == 'P'){
	 	phi = -1
	 	}
	 else{
	 	print("The Options Type is Wrong!")
	 	return
	 }
	if (vol == 0 or maturity == 0){
		return max(0.0, phi * (spot - strike))
		
	}
	if (strike == 0){
		print('The strike_price price cannot be zero')
		}
	d1 = (log(spot \strike) + (r + 0.5 * vol * vol) * maturity) \vol \sqrt(maturity)
	d2 = d1 - vol * sqrt(maturity)
    	return phi * spot * cdfNormal(0,1,phi * d1) - phi * strike * exp(-r * maturity) * cdfNormal(0,1,phi * d2)
}

def bs_vega(spot, strike, maturity, r, vol, q=0.0){
	/*
	 * spot：现货价，float:
	 * strike：行权价，float 
	 * maturity：剩余天数，float、int
	 * r：无风险利率，float
	 * vol：波动率，float
	*/
	d1 = (log(spot \ strike) + (r + vol * vol \ 2.) * maturity) \ (vol * sqrt(maturity))
	return spot * sqrt(maturity) * pdfNorm(d1)
}

def implied_volatility_Newton(options_price, options_type, spot, strike, maturity, r){
	/*
	 * 牛顿法实现期权隐含波动率计算
	 *  *options_price,期权价格，float:
	 * spot：现货价，float:
	 * strike：行权价，float 
	 * maturity：剩余天数，float、int
	 * r：无风险利率，float
	 * vol：波动率，float
	 * options_type:'期权类型 ，C' or 'P'
	 * 
	*/
	max_iterations= 100
	precision= 1.0e-5
	sigma = 0.5
	for (i in 1..max_iterations){
	    price = bs_price(options_type, spot, strike, maturity, r, sigma)
	    vega = bs_vega( spot, strike, maturity, r, sigma)
	    diff = options_price - price
	    if (abs(diff) < precision){
	        return sigma
	        }
	    sigma = sigma + diff / vega 
	}
	return sigma
}

def implied_volatility_Binary(options_price, spot, strike, r, maturity, options_type){
	/*
	 * 二分法实现期权隐含波动率计算
	 * options_price,期权价格，float:
	 * spot：现货价，float:
	 * strike：行权价，float 
	 * maturity：剩余天数，float、int
	 * r：无风险利率，float
	 * vol：波动率，float
	 * options_type:'期权类型 ，C' or 'P'
	*/
	sigma_min = 0.00001
	sigma_max = 1.000
	sigma_mid = (sigma_min + sigma_max) / 2
	sigma=each(bs_price{options_type, spot, strike, maturity, r, },[sigma_min,sigma_max,sigma_mid])
	if (options_price< sigma[0] || options_price > sigma[1]){
	    print('error, the price of option is beyond the limit')
	}
	diff = options_price - sigma[2]
	do{
		sigma=each(bs_price{options_type, spot, strike, maturity, r, },[sigma_mid, (sigma_min + sigma_max) / 2])
		diff = options_price - sigma[0]
		sigma_mid= (sigma_min + sigma_max) \2
		if (options_price > sigma[1]){
		    sigma_min =sigma_mid
		}
		else{
		    sigma_max = sigma_mid
		}
	}while (abs(diff) > 1e-6)    
	return  sigma_mid
}