DolphinDB 流计算在商品期货交易的应用：波动率计算与拟合

如前文所述，DolphinDB 提供时间邻近度关联左右表引擎，可以直接通过创建引擎的方式实现根据期权报价时间找最近的期货价格，如下所示：

ajEngine = createAsofJoinEngine(
  name="asofJoin",
  leftTable=ctpOptions,
  rightTable=ctpFutures,
  outputTable=futOptCrt,
  metrics=<[ctpOptions.InstrumentID, ctpFutures.assetCode, (ctpFutures.BidPrice1+ctpFutures.AskPrice1)\2, (ctpOptions.BidPrice1+ctpOptions.AskPrice1)\2, KPrice, CPMode, endDate]>,
  matchingColumn=[[`futureCode],[`InstrumentID]],
  timeColumn=`TradingTime,
  useSystemTime=false)
// subscribe topic
subscribeTable(tableName="ctpOptions", actionName="appendLeftStream", handler=getLeftStream(ajEngine), msgAsTable=true, offset=-1, hash=0)
subscribeTable(tableName="ctpFutures", actionName="appendRightStream", handler=getRightStream(ajEngine), msgAsTable=true, offset=-1, hash=1)

其中，期权行情为左表，期货行情为右表。对于每一笔期权行情数据，会从期货行情中获取临近的最新期货行情数据，并输出。读者可以通过参数delayedTime设置超时强制触发规则，详情请见引擎详细说明asof join 引擎。

brent 函数要求提供一个函数参数f，在一个给定区间 [a, b] 内求出函数 f(x) 的一个根x0，使得 f(x0)=0 ，所以在代码中除了实现 Black-76 的模型外，还定义了一个函数 loss_function，返回 Black-76 模型的定价与实际行情价格的差值，然后使用 brent 函数进行求解。计算代码如下：

/* 
 * @Brief: 计算隐波，black_76 公式
 * @Param: F: 标的期货的当前价格
 * @Param: K: 执行价
 * @Param: T: 剩余期限
 * @Param: r: 无风险利率
 * @Param: sigma: 波动率
 * @Param: option_type: 期权类型，1或-1，代表'C' or 'P'
*/
def black_76(F, K, T, r, sigma, option_type){
    d1 = (log(F \ K) + (0.5 * pow(sigma , 2) ) * T) \ (sigma * sqrt(T))
    d2 = d1 - sigma * sqrt(T)
    if (option_type == 1){
        return exp(-r * T) * (F * cdfNormal(0,1,d1) - K * cdfNormal(0,1,d2))
    }else if( option_type == -1){
        return exp(-r * T) * (K * cdfNormal(0,1,-d2) - F * cdfNormal(0,1,-d1))
    }else{
        return ("Option type must be 1 or -1.")
    }
}
/* 
 * @Brief: 计算隐波，Brentq 的求解函数
 * @Param: sigma: 波动率
 * @Param: F: 标的期货的当前价格
 * @Param: K: 执行价
 * @Param: T: 剩余期限
 * @Param: r: 无风险利率
 * @Param: option_type: 期权类型，1或-1，代表'C' or 'P'
 * @Param: market_price: 期权价格
*/
def loss_function(sigma, F, K, T, r, option_type, market_price){
    return black_76(F, K, T, r, sigma, option_type) - market_price
}
/**
 * @Brief: 计算隐波，Brentq 的求解函数
 * @Param: F: 标的期货的当前价格
 * @Param: K: 执行价
 * @Param: T: 剩余期限
 * @Param: r: 无风险利率
 * @Param: option_type: 期权类型 ，1或-1，代表'C' or 'P'
 * @Param: market_price: 期权价格
 * @Return: 
 * @Sample usage:
 */
def implied_volatility_Black76Brentq(F, K, T, r, option_type, market_price){
	iv = double(NULL)
	try{ 
			iv = brentq(f=loss_function, a=1e-6, b=2, xtol=2e-12, rtol=1e-9, maxIter=100, funcDataParam=[F, K, T, r, option_type, market_price])[1]
	}catch(ex){
			writeLog("calculate error: " + ex)
	}
	return iv
}

DolphinDB 的横截面引擎保存了所有分组的最新记录，有效解决部分合约因低频交易而缺少行情数据的问题。同时横截面引擎可以使用数据时间按照固定时间间隔触发。这里我们创建横截面引擎每隔数据时间1分钟，进行一次曲线构建与平滑。代码如下：

crossEngine = createCrossSectionalEngine(
  name="crossEngine",
  metrics=<[last(assetCode), buildIVCurve(trigger, futureCode, assetCode, futurePrice, KPrice, CPMode, iv)]>,
  dummyTable=ivResult,
  outputTable=ivCurve,
  keyColumn=`futureCode`KPrice`CPMode,
  triggeringPattern='dataInterval',
  triggeringInterval=60*1000,
  timeColumn=`trigger,
  useSystemTime=false,
  contextByColumn=`futureCode,
  outputElapsedMicroseconds=true)
filterFn = def (msg) { // 过滤掉 iv 无效的数据
	return (select * from msg where iv > 0)
}
subscribeTable(tableName="ivResult", actionName="buildCurve", offset=-1, handler=append!{crossEngine}, msgAsTable=true, batchSize=10000, throttle=1, hash=3, filter=filterFn);

通常我们使用流动性较好的虚值期权来构建波动率曲线。使用每分钟最新的期货价格来作为基准查找看涨看跌的虚值期权，代码如下：

strikePrices = KPrice[KPrice< lastPrice and CPMode == -1].join(KPrice[KPrice>= lastPrice and CPMode == 1])
strikeIvs = iv[KPrice< lastPrice and CPMode == -1].join(iv[KPrice>= lastPrice and CPMode == 1])
indexes = isort(strikePrices) // 排序
strikePrices = strikePrices[indexes]
strikeIvs = strikeIvs[indexes]

本文使用三次样条插值，先根据选取的所有虚值期权和隐含波动率构建模型，然后将中间的执行价格等分为50份，使用模型对这些价格进行插值，最终将原始的执行价及隐含波动率、平滑的执行价格及插值的隐含波动率，以 arrayVector 保存，代码如下：

n = 50 // 样本数量
if(strikePrices.size() >= 3 ){ // 波动率数据太少，就不做拟合
  // 构建微笑曲线执行价坐标
  curvePrices = linspace(min(strikePrices), max(strikePrices), n, true)[1]$int
  // 构建三次样条插值模型
  cs = cubicSpline(strikePrices, strikeIvs)
  // 插值获得微笑曲线的波动率数据
  curveIvs = cubicSplinePredict(cs, curvePrices)
  // 返回插值结果
  tableInsert(objByName(`interResult), ( max(trigger), first(futureCode), first(assetCode), arrayVector([strikePrices.size()], strikePrices), arrayVector([strikeIvs.size()], strikeIvs), arrayVector([curvePrices.size()], curvePrices), arrayVector([curveIvs.size()], curveIvs) ))
}else{
  prices, ivs = NULL, NULL
  if(strikePrices.count()>0) prices, ivs = arrayVector([strikePrices.size()], strikePrices), arrayVector([strikeIvs.size()], strikeIvs)
  tableInsert(objByName(`interResult), ( max(trigger), first(futureCode), first(assetCode), prices, ivs, NULL, NULL ))
}