如何提升 ETF 期權隱含波動率和希臘值的計算速度?

語言: CN / TW / HK

期權的隱含波動率可以反應市場對未來的預期,通常使用牛頓法和二分法來計算。這兩種方法都需要頻繁迭代,且迭代次數不能確定,核心程式碼無法向量化,因此只能通過迴圈來逼近求解。這就導致在期權相關計算中,隱含波動率往往容易成為效能的瓶頸。

DolphinDB 的計算邏輯使用指令碼語言編寫,但底層呼叫的是 C++ 程式碼,存在指令碼解釋的過程。為了提高指令碼的執行效率,DolphinDB 從 1.01 版本開始支援即時編譯(JIT)功能,特別適合於無法使用向量化運算但又對執行速度有極高要求的場景。

本教程將基於客戶的實際需求,以二分法計算 ETF 期權的隱含波動率及希臘值為例,為大家示範如何使用 DolphinDB 的 JIT 功能給計算過程加速,並與 C++ 原生程式碼進行了計算效能對比測試,結果表明 DolphinDB 指令碼計算耗時為 C++ 原生程式碼的1.5倍。

1. 資料表結構

1.1 期權日頻資料表

欄位

欄位型別

含義

tradedate

DATE

交易日期

sym

SYMBOL

標的程式碼

codes

SYMBOL

期權合約程式碼

closeprice

DOUBLE

日收盤價格

etf

SYMBOL

期權合成價格的兩個合約程式碼

etfprice

DOUBLE

期權合成價格

字串欄位使用 SYMBOL 型別和 STRING 型別儲存的差異,參考: 資料型別 — DolphinDB 2.0 文件的字串部分內容。

期權日頻資料表在 DolphinDB 中儲存時,建議在時間維度按年分割槽即可,建立庫表的程式碼如下:

login("admin", "123456")
dbName = "dfs://optionPrice"
tbName = "optionPrice"
if(existsDatabase(dbName)){
	dropDatabase(dbName)
}
db = database(dbName, RANGE, date(datetimeAdd(2000.01M,0..50*12,'M')))
colNames = `tradedate`sym`codes`closeprice`etf`etfprice
colTypes = [DATE, SYMBOL, SYMBOL, DOUBLE, SYMBOL, DOUBLE]
schemaTable = table(1:0, colNames, colTypes)
db.createPartitionedTable(table=schemaTable, tableName=tbName, partitionColumns=`tradedate)

讀取 50 ETF 日頻資料,程式碼如下:

data = select * from loadTable("dfs://optionPrice", "optionPrice") where sym =`510050

data 是表變數,具體資料內容如下圖所示:

通過 panel 函式對 data 進行透視操作,將窄表資料展成矩陣,生成期權日收盤價矩陣,程式碼如下:

closPriceWideMatrix = panel(data.codes, data.tradeDate, data.closePrice)

closPriceWideMatrix 是矩陣變數,具體資料內容如下圖所示:

通過 panel 函式對 data 進行透視操作,將窄表資料展成矩陣,生成期權合成價格矩陣,程式碼如下:

etfPriceWideMatrix = panel(data.codes, data.tradeDate, data.etfprice)

etfPriceWideMatrix 是矩陣變數,具體資料內容如下圖所示:

1.2 期權資訊

欄位

欄位型別

含義

code

STRING

期權合約程式碼

name

STRING

期權合約名稱

exemode

INT

期權型別(認購認沽型別)

exeprice

DOUBLE

行權價

startdate

DATE

開始日期

lastdate

DATE

結束日期

sym

SYMBOL

標的程式碼

exeratio

DOUBLE

乘數

exeprice2

DOUBLE

分紅後的行權價

dividenddate

DATE

分紅日

tradecode

STRING

交易程式碼

期權資訊資料在 DolphinDB 中儲存時,建議在證券程式碼維度按值分割槽即可,建立庫表的程式碼如下:

login("admin", "123456")
dbName = "dfs://optionInfo"
tbName = "optionInfo"
if(existsDatabase(dbName)){
	dropDatabase(dbName)
}
db = database(dbName, VALUE, `510050`510300)
colNames = `code`name`exemode`exeprice`startdate`lastdate`sym`exeratio`exeprice2`dividenddate`tradecode
colTypes = [STRING, STRING, INT, DOUBLE, DATE, DATE, SYMBOL, DOUBLE, DOUBLE, DATE, STRING]
schemaTable = table(1:0, colNames, colTypes)
db.createPartitionedTable(table=schemaTable, tableName=tbName, partitionColumns=`sym)

讀取資料:

contractInfo = select * from loadTable("dfs://optionInfo", "optionInfo") where sym =`510050

contractInfo 是表變數,具體資料內容如下圖所示:

1.3 交易日曆

欄位

欄位型別

含義

tradedate

DATE

交易日期

交易日曆存放在單列 csv 檔案中,可以使用 DolphinDB 的 loadText 函式介面直接讀取:

//交易日曆csv檔案路徑
tradingDatesAbsoluteFilename = "/hdd/hdd9/tutorials/jitAccelerated/tradedate.csv"
startDate = 2015.02.01
endDate = 2022.03.01
//讀取csv檔案
allTradingDates = loadText(tradingDatesAbsoluteFilename)
//生成交易日向量
tradingDates = exec tradedate from allTradingDates where tradedate<endDate and tradedate >startDate

allTradingDates 是表變數,可以使用 exec 函式將表中的某一列轉成向量變數,tradingDates 是向量變數。

2. 計算函式程式碼開發

2.1 隱含波動率

DolphinDB 指令碼語言需要先解釋再執行,計算密集的程式碼如果不能向量化,在指令碼層面使用 while 和 for 迴圈以及條件分支,就會比較耗時。期權隱含波動率計算的步驟,由於使用了上下限值迴圈逼近的二分法,正是需要 JIT 加速的計算型別。以下就是二分法逼近隱含波動率的程式碼:

@jit
def calculateD1JIT(etfTodayPrice, KPrice, r, dayRatio, HLMean){
	skRatio = etfTodayPrice / KPrice
	denominator = HLMean * sqrt(dayRatio)
	result = (log(skRatio) + (r + 0.5 * pow(HLMean, 2)) * dayRatio) / denominator
	return result
}

@jit
def calculatePriceJIT(etfTodayPrice, KPrice , r , dayRatio , HLMean , CPMode){
	testResult = 0.0
	if (HLMean <= 0){
		testResult = CPMode * (etfTodayPrice - KPrice)
		if(testResult<0){
			return 0.0
		}
		return testResult
	}
	d1 = calculateD1JIT(etfTodayPrice, KPrice, r, dayRatio, HLMean)
	d2 = d1 - HLMean * sqrt(dayRatio)
	price = CPMode * (etfTodayPrice * cdfNormal(0, 1, CPMode * d1) - KPrice * cdfNormal(0, 1, CPMode * d2) * exp(-r * dayRatio))
	return price
}

@jit
def calculateImpvJIT(optionTodayClose, etfTodayPrice, KPrice, r, dayRatio, CPMode){
	v = 0.0	
	high = 2.0
	low = 0.0
	do{
		if ((high - low) <= 0.00001){
			break
		}
		HLMean = (high + low) / 2.0
		if (calculatePriceJIT(etfTodayPrice, KPrice, r, dayRatio, HLMean, CPMode) > optionTodayClose){
			high = HLMean
		}
		else{
			low = HLMean
		}
	}
	while(true)
	v = (high + low) / 2.0
	return v
}

def calculateImpv(optionTodayClose, etfTodayPrice, KPrice, r, dayRatio, CPMode){
	originalShape = optionTodayClose.shape()
	optionTodayClose_vec = optionTodayClose.reshape()
	etfTodayPrice_vec = etfTodayPrice.reshape()
	KPrice_vec = KPrice.reshape()
	dayRatio_vec = dayRatio.reshape()
	CPMode_vec = CPMode.reshape()
	impvTmp = each(calculateImpvJIT, optionTodayClose_vec, etfTodayPrice_vec, KPrice_vec, r, dayRatio_vec, CPMode_vec)	
	impv = impvTmp.reshape(originalShape)	
	return impv
}

calculateImpvJIT 是計算隱含波動的核心程式碼,其入參 optionTodayClose, etfTodayPrice, KPrice, r, dayRatio, CPMode 都是標量物件,其呼叫的 calculatePriceJIT 函式和 calculateD1JIT 函式都通過 @jit 裝飾器的方式封裝成JIT函式,以達到加速計算的目的。

calculateImpv 是計算隱含波動的最終呼叫函式,其入參 optionTodayClose, etfTodayPrice, KPrice, dayRatio, CPMode 都是矩陣物件,其主要作用是把輸入和輸出進行矩陣和向量的轉換,以適不同函式的入參和輸出。在後面 delta, gamma, vega, theta 計算時,也會用到這些矩陣入參,這裡以 2015年2月16日的 50ETF 為例進行展示。

  • optionTodayClose

  • etfTodayPrice

  • KPrice

  • dayRatio

  • CPMode

2.2 delta

delta 表示期權價格對標的資產價格的變動率,即標的資產價格每變動一個單位,期權價格產生的變化。

delta 的計算可以方便地實現向量化計算,所以不需要呼叫 JIT 功能,其程式碼如下:

def calculateD1(etfTodayPrice, KPrice, r, dayRatio, HLMean){
	skRatio = etfTodayPrice / KPrice
	denominator = HLMean * sqrt(dayRatio)
	result = (log(skRatio) + (r + 0.5 * pow(HLMean, 2)) * dayRatio) / denominator
	return result
}

def cdfNormalMatrix(mean, stdev, X){
	originalShape = X.shape()
	X_vec = X.reshape()
	result = cdfNormal(mean, stdev, X_vec)
	return result.reshape(originalShape)
}

def calculateDelta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode){
	delta = iif(
			impvMatrix <= 0,
			0,
			0.01*etfTodayPrice*CPMode*cdfNormalMatrix(0, 1, CPMode * calculateD1(etfTodayPrice, KPrice, r, dayRatio, impvMatrix))
		)
	return delta
}

calculateDelta 是計算 delta 的最終呼叫函式,其入參 etfTodayPrice, KPrice, dayRatio, impvMatrix, CPMode 都是矩陣物件。

2.3 gamma

gamma 表示 delta 對於標的資產價格的變動率,即標的資產價格每變動一個單位,delta 值產生的變化。

gamma 的計算可以方便地實現向量化計算,所以不需要呼叫 JIT 功能,其程式碼如下:

def normpdf(x){
	return exp(-pow(x, 2)/2.0)/sqrt(2*pi)
}

def calculateD1(etfTodayPrice, KPrice, r, dayRatio, HLMean){
	skRatio = etfTodayPrice / KPrice
	denominator = HLMean * sqrt(dayRatio)
	result = (log(skRatio) + (r + 0.5 * pow(HLMean, 2)) * dayRatio) / denominator
	return result
}

def calculateGamma(etfTodayPrice, KPrice, r, dayRatio, impvMatrix){
	gamma = iif(
			impvMatrix <= 0,
			0,
			(normpdf(calculateD1(etfTodayPrice,  KPrice, r, dayRatio, impvMatrix)) \ (etfTodayPrice * impvMatrix * sqrt(dayRatio))) * pow(etfTodayPrice, 2) * 0.0001
		)	
	return gamma
}

calculateGamma 是計算 gamma 的最終呼叫函式,其入參 etfTodayPrice, KPrice, dayRatio, impvMatrix 都是矩陣物件。

2.4 vega

vega 表示波動率單位變動對期權價格產生的變化。

vega 的計算可以方便地實現向量化計算,所以不需要呼叫 JIT 功能,其程式碼如下:

def normpdf(x){
	return exp(-pow(x, 2)/2.0)/sqrt(2*pi)
}

def calculateD1(etfTodayPrice, KPrice, r, dayRatio, HLMean){
	skRatio = etfTodayPrice / KPrice
	denominator = HLMean * sqrt(dayRatio)
	result = (log(skRatio) + (r + 0.5 * pow(HLMean, 2)) * dayRatio) / denominator
	return result
}

def calculateVega(etfTodayPrice, KPrice, r, dayRatio, impvMatrix){
	vega = iif(
			impvMatrix <= 0,
			0,
			etfTodayPrice * normpdf(calculateD1(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)) * sqrt(dayRatio)
		)
	return vega \ 100.0
}

calculateVega 是計算 vega 的最終呼叫函式,其入參 etfTodayPrice, KPrice, dayRatio, impvMatrix 都是矩陣物件。

2.5 theta

theta 表示時間流逝對期權價格產生的變化,即每減少一天,期權價格的變化值。

theta 的計算可以方便地實現向量化計算,所以不需要呼叫 JIT 功能,其程式碼如下:

def calculateD1(etfTodayPrice, KPrice, r, dayRatio, HLMean){
	skRatio = etfTodayPrice / KPrice
	denominator = HLMean * sqrt(dayRatio)
	result = (log(skRatio) + (r + 0.5 * pow(HLMean, 2)) * dayRatio) / denominator
	return result
}

def normpdf(x){
	return exp(-pow(x, 2)/2.0)/sqrt(2*pi)
}

def cdfNormalMatrix(mean, stdev, X){
	originalShape = X.shape()
	X_vec = X.reshape()
	result = cdfNormal(mean, stdev, X_vec)
	return result.reshape(originalShape)
}


def calculateTheta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode){
	annualDays = 365
	d1 = calculateD1(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)
	d2 = d1 - impvMatrix * sqrt(dayRatio)
	theta = (-etfTodayPrice * normpdf(d1) * impvMatrix \ (2 * sqrt(dayRatio)) - CPMode * r * KPrice * exp(-r * dayRatio) *cdfNormalMatrix(0, 1, CPMode * d2)) \ annualDays
	result = iif(impvMatrix<= 0, 0, theta)	
	return result
}

calculateTheta 是計算 theta 的最終呼叫函式,其入參 etfTodayPrice, KPrice, dayRatio, impvMatrix, CPMode 都是矩陣物件。

2.6 單日計算函式

開發完最核心的計算函式後,可以自定義一個單日計算函式,計算指定日期的隱含波動率和希臘值,其程式碼如下:

def calculateOneDayGreek(closPriceWideMatrix, etfPriceWideMatrix, contractInfo, targetDate){
	targetDate_vec = [targetDate]
	r = 0
	optionTodayClose = getTargetDayOptionClose(closPriceWideMatrix, targetDate_vec)
	validContractsToday = optionTodayClose.columnNames()
	etfTodayPrice = getTargetDayEtfPrice(etfPriceWideMatrix, targetDate_vec)
	KPrice, dayRatio, CPMode = getTargetDayContractInfo(contractInfo, validContractsToday, targetDate_vec)
	impvMatrix = calculateImpv(optionTodayClose, etfTodayPrice, KPrice, r, dayRatio, CPMode)
	deltaMatrix = calculateDelta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode)\(etfTodayPrice*0.01)
	gammaMatrix = calculateGamma(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)\(pow(etfTodayPrice, 2) * 0.0001)
	vegaMatrix = calculateVega(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)
	thetaMatrix = calculateTheta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode)
	todayTable = table(validContractsToday as optionID, impvMatrix.reshape() as impv, deltaMatrix.reshape() as delta, gammaMatrix.reshape() as gamma, vegaMatrix.reshape() as vega, thetaMatrix.reshape() as theta)
	todayTable["tradingDate"] = targetDate
	todayTable.reorderColumns!(["optionID", "tradingDate"])
	return todayTable
}

calculateOneDayGreek 的入參 closPriceWideMatrix, etfPriceWideMatrix 是矩陣物件(參考第1章節的讀取資料),contractInfo 是表物件(參考第1章節的讀取資料),targetDate 是標量物件。

calculateOneDayGreek 函式還呼叫了 getTargetDayOptionClose 函式, getTargetDayEtfPrice 函式和getTargetDayContractInfo 函式,呼叫的目的是從全量資料中獲取計算當日的有效資訊,程式碼如下:

/*
 * 按合約和交易日在期權日頻收盤價矩陣中尋找對應價格
 */
def getTargetDayOptionClose(closPriceWideMatrix, targetDate){
	colNum = closPriceWideMatrix.colNames().find(targetDate)
	return closPriceWideMatrix[colNum].transpose().dropna(byRow = false)
}

/*
 * 按合約和交易日在期權合成期貨價格矩陣中尋找對應價格
 */
def getTargetDayEtfPrice(etfPriceWideMatrix, targetDate){
	colNum = etfPriceWideMatrix.colNames().find(targetDate)
	return etfPriceWideMatrix[colNum].transpose().dropna(byRow = false)
}

/*
 * 根據合約和交易日在期權資訊表中尋找 KPrice, dayRatio, CPMode
 */
def getTargetDayContractInfo(contractInfo, validContractsToday, targetDate){
	targetContractInfo = select code, exemode, exeprice, lastdate, exeprice2, dividenddate, targetDate[0] as tradingDate from contractInfo where Code in validContractsToday
	KPrice = exec iif(tradingDate<dividenddate, exeprice2, exeprice) from targetContractInfo pivot by tradingDate, code
	dayRatio = exec (lastdate-tradingDate)\365.0 from targetContractInfo pivot by tradingDate, Code
	CPMode = exec exemode from targetContractInfo pivot by tradingDate, Code
	return KPrice, dayRatio, CPMode
}

calculateOneDayGreek 函式的具體使用方法會在下一章說明。

2.7 多日平行計算函式

自定義單日計算函式後,可以再自定義一個多日平行計算函式,其程式碼如下:

def calculateAll(closPriceWideMatrix, etfPriceWideMatrix, contractInfo, tradingDates, mutable result){
	calculator = partial(calculateOneDayGreek, closPriceWideMatrix, etfPriceWideMatrix, contractInfo)
	timer{
		allResult = ploop(calculator, tradingDates)
	}
	for(oneDayResult in allResult){
		append!(result, oneDayResult)
	}	
}

calculateAll 是自定義的多日平行計算函式,主要用到了 DolphinDB 內建的 partial 部分應用函式ploop 平行計算函式,直接傳入要並行的函式和入參,不必像其他語言一樣先定義執行緒池/程序池。calculateAll 函式的具體使用方法會在下一章說明。

3. 計算效能測試

測試的交易日範圍從2015年2月到2022年3月,實際交易日1729天。測試的期權品種是 50 ETF,證券程式碼為510050,涉及期權合約共3124個。

3.1 測試環境

  • CPU 型別:Intel(R) Xeon(R) Silver 4216 CPU @ 2.10GHz

  • 邏輯 CPU 總數:8

  • 記憶體:64GB

  • OS:64位 CentOS Linux 7 (Core)

3.2 單日計算效能測試

單執行緒單日計算效能測試程式碼如下:

//定義單日效能測試函式
def testOneDayPerformance(closPriceWideMatrix, etfPriceWideMatrix, contractInfo, targetDate){
	targetDate_vec = [targetDate]
	r = 0
	optionTodayClose = getTargetDayOptionClose(closPriceWideMatrix, targetDate_vec)
	validContractsToday = optionTodayClose.columnNames()
	etfTodayPrice = getTargetDayEtfPrice(etfPriceWideMatrix, targetDate_vec)
	KPrice, dayRatio, CPMode = getTargetDayContractInfo(contractInfo, validContractsToday, targetDate_vec)
	timer{
		impvMatrix = calculateImpv(optionTodayClose, etfTodayPrice, KPrice, r, dayRatio, CPMode)
		deltaMatrix = calculateDelta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode)\(etfTodayPrice*0.01)
		gammaMatrix = calculateGamma(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)\(pow(etfTodayPrice, 2) * 0.0001)
		vegaMatrix = calculateVega(etfTodayPrice, KPrice, r, dayRatio, impvMatrix)
		thetaMatrix = calculateTheta(etfTodayPrice, KPrice, r, dayRatio, impvMatrix, CPMode)
	}
	todayTable = table(validContractsToday as optionID, impvMatrix.reshape() as impv, deltaMatrix.reshape() as delta, gammaMatrix.reshape() as gamma, vegaMatrix.reshape() as vega, thetaMatrix.reshape() as theta)
	todayTable["tradingDate"] = targetDate
	todayTable.reorderColumns!(["optionID", "tradingDate"])
	return todayTable
}
//執行單日效能測試函式
oneDay = testOneDayPerformance(closPriceWideMatrix, etfPriceWideMatrix, contractInfo, 2022.02.28)

測試結果如下:

  • 計算日期為 2022年2月28日

  • 測試的期權品種是 50 ETF,涉及期權合約共136個

  • DolphinDB 指令碼計算總耗時為2.1 ms

  • C++ 原生程式碼計算總耗時為1.02 ms

3.3 多日平行計算效能測試

多日平行計算效能測試程式碼如下:

//建立儲存計算結果的表變數
result = table(
		array(SYMBOL, 0) as optionID,
		array(DATE, 0) as tradingDate,
		array(DOUBLE, 0) as impv,
		array(DOUBLE, 0) as delta,
		array(DOUBLE, 0) as gamma,
		array(DOUBLE, 0) as vega,
		array(DOUBLE, 0) as theta
	)
//執行多日平行計算函式
calculateAll(closPriceWideMatrix, etfPriceWideMatrix, contractInfo, tradingDates, result)

測試結果如下:

  • 計算日期為 2015年2月到2022年3月

  • 測試的期權品種是 50 ETF,涉及期權合約共3124個

  • 計算的並行度為8,測試環境的8個 CPU 滿負荷執行

  • DolphinDB 指令碼計算總耗時為300 ms

  • C++ 原生程式碼計算總耗時為200 ms

計算過程中的 CPU 使用率:

4. 總結

本教程中期權隱含波動率的計算使用了 JIT 功能提速,其餘希臘值的計算使用了向量化計算。我們測試了2015年2月到2022年3月 50 ETF 所有期權合約的隱含波動率和希臘值的計算效能,在8個 CPU 滿負荷執行下,DolphinDB 指令碼計算總耗時為300 ms,C++ 原生程式碼計算總耗時為200 ms,耗時相差50%左右。

關於 DolphinDB JIT 的更多詳細特性,可以參考 DolphinDB JIT教程

通過 DolphinDB 下載連線下載 DolphinDB server 進行測試的時候,必須選擇包含 JIT 功能的安裝包,如下圖所示:

附錄

呼叫計算函式的計算指令碼.txt

期權隱含波動率和希臘值計算函式.txt