# 通过 DolphinDB JIT 功能加速计算 ETF 期权隐含波动率和希腊值 期权的隐含波动率可以反应市场对未来的预期,通常使用牛顿法和二分法来计算。这两种方法都需要频繁迭代,且迭代次数不能确定,核心代码无法向量化,因此只能通过循环来逼近求解。这就导致在期权相关计算中,隐含波动率往往容易成为性能的瓶颈。 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 文档](https://www.dolphindb.cn/cn/help/DataTypesandStructures/DataTypes/index.html)的字符串部分内容。 期权日频数据表在 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 函数](https://www.dolphindb.cn/cn/help/FunctionsandCommands/FunctionReferences/p/panel.html)对 data 进行透视操作,将窄表数据展成矩阵,生成期权日收盘价矩阵,代码如下: ``` closPriceWideMatrix = panel(data.codes, data.tradeDate, data.closePrice) ``` closPriceWideMatrix 是矩阵变量,具体数据内容如下图所示: 通过 [panel 函数](https://www.dolphindb.cn/cn/help/FunctionsandCommands/FunctionReferences/p/panel.html)对 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 函数](https://www.dolphindb.cn/cn/help/FunctionsandCommands/FunctionReferences/l/loadText.html)接口直接读取: ``` //交易日历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 tradedatestartDate ``` 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 ### 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教程](https://gitee.com/dolphindb/Tutorials_CN/blob/master/jit.md)。 通过 [DolphinDB 下载连接](https://www.dolphindb.cn/alone/alone.php?id=10)下载 DolphinDB server 进行测试的时候,必须选择包含 JIT 功能的安装包,如下图所示: ## 附录 [调用计算函数的计算脚本.txt](./script/IV_Greeks_Calculation_for_ETF_Options_Using_JIT/calculation_scripts.txt) [期权隐含波动率和希腊值计算函数.txt](./script/IV_Greeks_Calculation_for_ETF_Options_Using_JIT/Option_IV_and_GreekS_Calculation_Functions.txt)