Feature Engineering for Stock Volatility Prediction

The following indicators are calculated based on the 10 highest bid prices and 10 lowest ask prices mentioned above.

• wap: weighted average price

  
  
  

• priceSpread: bid-ask spread

  
  
  

• bidSpread: difference between the lowest and second lowest bid prices

  
  
  

• offerSpread: difference between the highest and second highest ask prices

  
  
  

• totalVolume: total size of the 10 highest bids and 10 lowest asks

  
  
  

• volumeImbalance: absolute value of the difference between the sizes of 10 highest bid prices and 10 lowest ask prices

  
  
  

• logReturnOffer: log return of ask prices

  
  
  

• logReturnBid: log return of bid prices

  
  
  

• wapBalance: absolute value of wap difference

  
  
  

• logReturnWap: log return of wap

  
  
  

The derived features are generated by aggregating over the indicators with 10-minute windows.

The following aggregate functions are used:

We will predict the realized volatility:

Considering the features are time-sensitive, we form 4 overlapping windows based on the 10-minute window (0-600s, 150-600s, 300-600s, and 450-600s) and conduct calculations on all of them.

Finally, 676 dimensions of derived features are generated.

The following part in function featureEngineering calculates the indicators:

 wap = (BidPrice * OfferOrderQty + BidOrderQty * OfferPrice) \ (BidOrderQty + OfferOrderQty)
 wapBalance = abs(wap[0] - wap[1])
 priceSpread = (OfferPrice[0] - BidPrice[0]) \ ((OfferPrice[0] + BidPrice[0]) \ 2)
 BidSpread = BidPrice[0] - BidPrice[1]
 OfferSpread = OfferPrice[0] - OfferPrice[1]
 totalVolume = OfferOrderQty.rowSum() + BidOrderQty.rowSum()
 volumeImbalance = abs(OfferOrderQty.rowSum() - BidOrderQty.rowSum())
 LogReturnWap = logReturn(wap)
 LogReturnOffer = logReturn(OfferPrice)
 LogReturnBid = logReturn(BidPrice)

In Python, we can pass a dictionary (where key is column names and value is a list of aggregate functions) to the pandas function groupby.agg to implement aggregation on selected columns.

In DolphinDB, we use a user-defined function createAggMetaCode to convert the above dictionary to SQL statements with metaprogramming.

 def createAggMetaCode(aggDict){
     metaCode = []
     metaCodeColName = []
     for(colName in aggDict.keys()){
         for(funcName in aggDict[colName])
         {
             metaCode.append!(sqlCol(colName, funcByName(funcName), colName + `_ + funcName$STRING))
             metaCodeColName.append!(colName + `_ + funcName$STRING)
         }
     }
     return metaCode, metaCodeColName$STRING
 }
 ​
 features = {
     "DateTime":[`count]
 }
 for( i in 0..9)
 {
     features["Wap"+i] = [`sum, `mean, `std]
     features["LogReturn"+i] = [`sum, `realizedVolatility, `mean, `std]
     features["LogReturnOffer"+i] = [`sum, `realizedVolatility, `mean, `std]
     features["LogReturnBid"+i] = [`sum, `realizedVolatility, `mean, `std]
 }
 features["WapBalance"] = [`sum, `mean, `std]
 features["PriceSpread"] = [`sum, `mean, `std]
 features["BidSpread"] = [`sum, `mean, `std]
 features["OfferSpread"] = [`sum, `mean, `std]
 features["TotalVolume"] = [`sum, `mean, `std]
 features["VolumeImbalance"] = [`sum, `mean, `std]
 aggMetaCode, metaCodeColName = createAggMetaCode(features)

It returns a vector of metacode and the corresponding column names:

The function featureEngineering generates a table with indicators and renamed columns, which facilitates the calculation of derived features with metaprogramming.

subTable = table(DateTime as `DateTime, BidPrice, BidOrderQty, OfferPrice, OfferOrderQty, wap, wapBalance, priceSpread, BidSpread, OfferSpread, totalVolume, volumeImbalance, LogReturnWap, LogReturnOffer, LogReturnBid)
colNum = 0..9$STRING
colName = `DateTime <- (`BidPrice + colNum) <- (`BidOrderQty + colNum) <- (`OfferPrice + colNum) <- (`OfferOrderQty + colNum) <- (`Wap + colNum) <- `WapBalance`PriceSpread`BidSpread`OfferSpread`TotalVolume`VolumeImbalance <- (`LogReturn + colNum) <- (`LogReturnOffer + colNum) <- (`LogReturnBid + colNum)
subTable.rename!(colName) 

Note: The operator "<-" merges multiple fields into a vector.

Finally, the metacode is passed to function featureEngineering and the derived features are calculated based on the table generated in the last step. It returns a matrix with 676 columns.

 subTable['BarDateTime'] = bar(subTable['DateTime'], 10m)
 result = sql(select = aggMetaCode, from = subTable).eval().matrix()
 result150 = sql(select = aggMetaCode, from = subTable, where = <time(DateTime) >= (time(BarDateTime) + 150*1000) >).eval().matrix()
 result300 = sql(select = aggMetaCode, from = subTable, where = <time(DateTime) >= (time(BarDateTime) + 300*1000) >).eval().matrix()
 result450 = sql(select = aggMetaCode, from = subTable, where = <time(DateTime) >= (time(BarDateTime) + 450*1000) >).eval().matrix()
 return concatMatrix([result, result150, result300, result450])

Part of the output table:

See the full scripts of calculating 676 derived features with DolphinDB metaprogramming:

• Batch processing in DolphinDB: storing 10 highest bids and 10 lowest asks with multiple columns
• Batch processing in DolphinDB (with array vectors): storing 10 levels of quotes with array vectors

The following script first excludes records with NULL values from the table. The metric LogReturn0_realizedVolatility is used to label the results. The table is adjusted to be passed to the XGBoost method.

//delete records with NULL values
 result = result[each(isValid, result.values()).rowAnd()]
 ​
 result_input = copy(result)
 ​
//Use LogReturn0_realizedVolatility as the label
 label = result[`LogReturn0_realizedVolatility]
 ​
//convert the data type of column SecurityID to INT
 result_input.update!(`SecurityID_int, int(result[`SecurityID]))
 ​
//remove columns that are not needed
 result_input.dropColumns!(`SecurityID`DateTime`LogReturn0_realizedVolatility)

Note: This project aims to predict the realized volatility in the next 10 minutes. As Wap_0 is the closest to the stock price, we use LogReturn0_realizedVolatility (which is calculated based on Wap_0) as the label.

We do not set validation dataset in this project. The data is split into training set and testing set with a 70-30 split.

 def trainTestSplit(x, testRatio) {
     xSize = x.size()
     testSize =( xSize * (1-testRatio))$INT
     return x[0: testSize], x[testSize:xSize]
 }
 ​
 Train_x, Test_x = trainTestSplit(result_input, 0.3)
 Train_y, Test_y = trainTestSplit(label, 0.3)

DolphinDB XGBoost plugin provides the following 4 methods:

• xgboost::train(): train a model
• xgboost::predict: make predictions with a model
• xgboost::saveModel: save a model to disk
• xgboost::loadModel: load a model from disk

The following script is used for model training and evaluation in DolphinDB:

 //define the validation metrics RMSPE
 def RMSPE(a,b)
 {
     return sqrt( sum( ((a-b)\a)*((a-b)\a) ) \a.size()  )
 }
 ​
 //define parameters for XGBoost training
 params = {
     objective: 'reg:squarederror',
     colsample_bytree: 0.8,
     subsample: 0.8,
     min_child_weight: 1,
     max_leaves:128,
     eta: 0.1,
     max_depth:10,
     eval_metric : 'rmse'
     }
 ​
 //model training
 model_1 = xgboost::train(Train_y ,Train_x, params, 500)
 ​
 //predict with the testing set and calculate RMPSE
 y_pred = xgboost::predict(model_1, Test_x)
 print('RMSPE='+RMSPE(Test_y, y_pred))

Result:

 RMSPE: 0.559
 Training Time: 1m 3s 327ms

Note: The params defined above is manually adjusted.

Save and load the model:

//save the model to modelSavePath, you can modify the path as you need
xgboost::saveModel(model_1, modelSavePath)

//load the model from modelSavePath, you can modify the path as you need
model = xgboost::loadModel(modelSavePath)

Prediction Performance Evaluation

The above diagram depicts the streaming process:

1. Ingest real-time stream data into table snapshotStream via DolphinDB API;
2. Subscribe to the table and ingest the data into DolphinDB time-series stream engine where calculations are conducted in a sliding window of 10 minutes (with step of 10 minutes).

The following is the main script for stream processing:

• Create the following stream tables:

  • snapshotStream: to store the snapshot data

  • aggrFeatures10min: the output table of feature engineering

  • result10min: the output table of predictions

    name = `SecurityID`DateTime`PreClosePx`OpenPx`HighPx`LowPx`LastPx`TotalVolumeTrade`TotalValueTrade`BidPrice0`BidPrice1`BidPrice2`BidPrice3`BidPrice4`BidPrice5`BidPrice6`BidPrice7`BidPrice8`BidPrice9`BidOrderQty0`BidOrderQty1`BidOrderQty2`BidOrderQty3`BidOrderQty4`BidOrderQty5`BidOrderQty6`BidOrderQty7`BidOrderQty8`BidOrderQty9`OfferPrice0`OfferPrice1`OfferPrice2`OfferPrice3`OfferPrice4`OfferPrice5`OfferPrice6`OfferPrice7`OfferPrice8`OfferPrice9`OfferOrderQty0`OfferOrderQty1`OfferOrderQty2`OfferOrderQty3`OfferOrderQty4`OfferOrderQty5`OfferOrderQty6`OfferOrderQty7`OfferOrderQty8`OfferOrderQty9
    
    type =`SYMBOL`TIMESTAMP`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`INT`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`INT`INT`INT`INT`INT`INT`INT`INT`INT`INT`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`DOUBLE`INT`INT`INT`INT`INT`INT`INT`INT`INT`INT
    
    share streamTable(100000:0, name, type) as snapshotStream
    share streamTable(100000:0 , `DateTime`SecurityID <- metaCodeColName <- (metaCodeColName+"_150") <- (metaCodeColName+"_300") <- (metaCodeColName+"_450"),`TIMESTAMP`SYMBOL <- take(`DOUBLE, 676)) as aggrFeatures10min
    share streamTable(100000:0 , `Predicted`SecurityID`DateTime, `FLOAT`SYMBOL`TIMESTAMP) as result10min

• Register a time-series stream engine

  metrics=sqlColAlias(<featureEngineering(DateTime,
  		matrix(BidPrice0,BidPrice1,BidPrice2,BidPrice3,BidPrice4,BidPrice5,BidPrice6,BidPrice7,BidPrice8,BidPrice9),
  		matrix(BidOrderQty0,BidOrderQty1,BidOrderQty2,BidOrderQty3,BidOrderQty4,BidOrderQty5,BidOrderQty6,BidOrderQty7,BidOrderQty8,BidOrderQty9),
  		matrix(OfferPrice0,OfferPrice1,OfferPrice2,OfferPrice3,OfferPrice4,OfferPrice5,OfferPrice6,OfferPrice7,OfferPrice8,OfferPrice9),
  		matrix(OfferOrderQty0,OfferOrderQty1,OfferOrderQty2,OfferOrderQty3,OfferOrderQty4,OfferOrderQty5,OfferOrderQty6,OfferOrderQty7,OfferOrderQty8,OfferOrderQty9), aggMetaCode)>, metaCodeColName <- (metaCodeColName+"_150") <- (metaCodeColName+"_300") <- (metaCodeColName+"_450"))
  
  createTimeSeriesEngine(name="aggrFeatures10min", windowSize=600000, step=600000, metrics=metrics, dummyTable=snapshotStream, outputTable=aggrFeatures10min, timeColumn=`DateTime, useWindowStartTime=true, keyColumn=`SecurityID)

• First subscription: subscribe to table snapshotStream

  subscribeTable(tableName="snapshotStream", actionName="aggrFeatures10min", offset=-1, handler=getStreamEngine("aggrFeatures10min"), msgAsTable=true, batchSize=2000, throttle=1, hash=0, reconnect=true)

• Second subscription: subscribe to table aggrFeatures10min

  def predictRV(mutable result10min, model, mutable msg){
  	startTime = now()
  	temp_table = select SecurityID, DateTime from msg
  	msg.update!(`SecurityID_int, int(msg[`SecurityID])).dropColumns!(`SecurityID`DateTime`LogReturn0_realizedVolatility)
  	Predicted = xgboost::predict(model , msg)
  	temp_table_2 = table(Predicted, temp_table)
  	result10min.append!(temp_table_2)
  }
  subscribeTable(tableName="aggrFeatures10min", actionName="predictRV", offset=-1, handler=predictRV{result10min, model}, msgAsTable=true, hash=1, reconnect=true)

  As DolphinDB has implemented a unified framework of batch and stream processing, the function featureEngineering called in the metrics defined above is identical to the user-defined function featureEngineering used in batch processing.

  After constructing the streaming framework, the test data is replayed into a shared stream table. The predicted realized volatility is output as follows:

  
  
  

In this section, we measure the computational latency in the time-series engine. It consists of 2 parts: the execution time of calculating aggregate features and of predicting realized volatility.

• Measure the execution time of calculating aggregate features:

  timer getStreamEngine('aggrFeatures10min').append!(data)

• Measure the execution time of predicting realized volatility:

  test_x  = select * from aggrFeatures10min
  timer{
  temp_table = select SecurityID, DateTime from test_x
  test_x.update!(`SecurityID_int, int(test_x[`SecurityID])).dropColumns!(`SecurityID`DateTime`LogReturn0_realizedVolatility)
  Predicted = xgboost::predict(model , test_x)
  }

• Results

  Number of Securities	Number of Records (in 10 min)	Calculating Aggregate Features (ms)	Predicting Realized Volatility (ms)	Total Latency (ms)
  1	201	22	3	25
  10	2011	92	3	95
  20	4020	162	4	168
  30	6030	257	5	262
  40	8040	321	6	327
  50	10054	386	7	393

See the full script of stream processing in DolphinDB:

Stream processing in DolphinDB

Test data (3-second snapshots of level 2 tick data)