如何提高在 DolphinDB 中计算希腊语的性能?



我正在使用DolphinDB来计算希腊语,我以矢量化的方式编写了它,性能相当不错。但我无法以矢量化的方式实现隐含波动率,这使得性能非常差。如何提高以下实现的性能?

def GBlackScholes(future_price, strike, input_ttm, risk_rate, b_rate, input_vol, is_call) {
ttm = input_ttm + 0.000000000000001;
vol = input_vol + 0.000000000000001;
d1 = (log(future_price/strike) + (b_rate + vol*vol/2) * ttm) / (vol * sqrt(ttm));
d2 = d1 - vol * sqrt(ttm);
if (is_call) {
return future_price * exp((b_rate - risk_rate) * ttm) * cdfNormal(0, 1, d1) - strike * exp(-risk_rate*ttm) * cdfNormal(0, 1, d2);
} else {
return strike * exp(-risk_rate*ttm) * cdfNormal(0, 1, -d2) - future_price * exp((b_rate - risk_rate) * ttm) * cdfNormal(0, 1, -d1);
}
}
def ImpliedVolatility(future_price, strike, ttm, risk_rate, b_rate, option_price, is_call) {
high=5.0;
low = 0.0;
do {
if (GBlackScholes(future_price, strike, ttm, risk_rate, b_rate, (high+low)/2, is_call) > option_price) {
high = (high+low)/2;
} else {
low = (high + low) /2;
}
} while ((high-low) > 0.00001);
return (high + low) /2;
}

如何提高以下实现的性能?

矢量化处理 ?

这部分有点神秘 - 原谅一个老量化没有很好地阅读这一点 - 关于哪些参数不打算成为标量的信息不存在,因此分析是基于明确呈现的信息片段。

免责声明
虽然我知道 DolphinDB 没有发布一个三元运算符(...)?(...):(...)可以在公共 API 中使用,但我觉得下面表达的想法是清晰而合理的。

战斗性能 ?

如果认真谈论性能,首先让我们回顾一下,并主要避免在上面提出的代码中发生的任何重复的重新计算:

def GBlackScholes( future_price,
strike,
input_ttm,
risk_rate,
b_rate,
input_vol,               // <------------[VAR]:: <-- ( high+low )/2
is_call                 //
) {                    //
ttm = input_ttm + 0.000000000000001; //-do-while-(CONST)
vol = input_vol + 0.000000000000001;//--do-while--------[VAR]
//
d1  = ( log( future_price )       //----do-while-(CONST)
- log( strike )            //-----do-while-(CONST)
+ b_rate                  //------do-while-(CONST)
+ vol*vol/2              //-------do-while--------[VAR]
* ttm                   //--------do-while-(CONST)
) / ( vol            //---------do-while--------[VAR]
* sqrt( ttm )   //----------do-while-(CONST)
);           //
d2  = ( d1                //------------do-while--------[VAR]
- vol              //-------------do-while--------[VAR]
* sqrt( ttm )     //--------------do-while-(CONST)
);             //              ++---------------[VAR]
//               ||     .________________________________________________.
// -----------[VAR]-?-( cdfNormal(--------vv-) * [                                                ]--do-while-(CONST)
return ( is_call ? ( cdfNormal( 0, 1,  d1 ) * future_price * exp( ( b_rate - risk_rate ) * ttm )
- cdfNormal( 0, 1,  d2 ) * strike       * exp(           -risk_rate   * ttm )
)
: ( cdfNormal( 0, 1, -d2 ) * strike       * exp(           -risk_rate   * ttm )
- cdfNormal( 0, 1, -d1 ) * future_price * exp( ( b_rate - risk_rate ) * ttm )
)
);
}

一个稍微好一点的GBlackScholes_WHILEd()函数 - 这为每个do{}while节省了~ 22xfloat-OP(其中一些非常昂贵(:

def GBlackScholes_WHILEd( vol_,         // <--------------------------[VAR]:: ( high+low )/2 + 0.000000000000001;
V1,           // <--------------------------[VAR]:: vol_ * C3
is_call,
ttm_,         // do-while-(CONST)
C1, C2, C3,   // do-while-(CONST)
R1, R2        // do-while-(CONST)
) {
d1  = ( C1
+ C2
* vol_*vol_ //--------------------do-while--------[VAR]
)
/ V1;       //--------------------do-while--------[VAR]
d2  = ( d1 //---------------------------do-while--------[VAR]
- V1 //---------------------------do-while--------[VAR]
);
//    --[VAR]--- ? ( cdfNormal(------[VAR]) * <________________________________________________>--do-while-(CONST)
return ( is_call ? ( cdfNormal( 0, 1,  d1 ) * R1
- cdfNormal( 0, 1,  d2 ) * R2
)
: ( cdfNormal( 0, 1, -d2 ) * R2
- cdfNormal( 0, 1, -d1 ) * R1
)
);
}

最后:

最有效的ImpliedVolatility()函数甚至完全避免了每个循环的所有调用签名处理,做了一些代数并保持在可实现的性能优势:

def ImpliedVolatility( future_price,
strike,
ttm,
risk_rate,
b_rate,
option_price,
is_call
) {
high = 5.0;                                                 // IS THIS A UNIVERSALLY SAFE & TRUE SUPREME - i.e. SAFELY ABOVE ALL POSSIBLE OPTIONS ?
low  = 0.0;
ttm_ = ttm + 0.000000000000001;                             // do-while-(CONST)     1x fADD
C1   = log(future_price ) - log( strike ) + b_rate;         // do-while-(CONST)     1x fADD  1x fDIV  1x fLOG  1x fNEG
C2   =     ( ttm_ ) / 2;                                    // do-while-(CONST)              1x fDIV
C3   = sqrt( ttm_ );                                        // do-while-(CONST)                       1x fSQRT
R1   = future_price * exp( ( b_rate - risk_rate ) * ttm_ ); // do-while-(CONST)     1x fADD  2x fMUL  1x fEXP  1x fNEG
R2   = strike       * exp(           -risk_rate   * ttm_ ); // do-while-(CONST)              2x fMUL  1x fEXP  1x fNEG
U4   = C2   - ttm_;                                         // do-while-(CONST)     1x fADD                    1x fNEG
U5   = ttm_ - C2;                                           // do-while-(CONST)     1x fADD                    1x fNEG
U3inv= 1./ C3;                                              // do-while-(CONST)              1x fDIV
// ------------------------------------------------------------// -----------------------------------------------------------------------------------------------
if ( is_call ) {                                            // do-while-RE-TESTING: AVOIDED REPETITIVE per-loop COSTS of TESTING THE VERY THE SAME
// -----------------------------------------------------------------------------------------------
do {
mid   = ( high + low ) / 2;       // cheapest  do-while per-loop-[VAR]-update
vol_  = mid + 0.000000000000001; //  cheapest  do-while per-loop-[VAR]-update
vol_2 = vol_ * vol_;            //   cheapest  do-while per-loop-[VAR]-update
/* ---------------------------------------------------------------------------------------------------------------------------------------------
HAS EVOLVED FROM THE ORIGINAL FORMULATION + AVOIDED REPETITIVE per-loop COSTS of 20+ expensive float OPs fully wasted,all in do-while-(CONST)
+ AVOIDED REPETITIVE per-loop COSTS of all the CALL fun() STACK MANIPULATIONS AND RELATED OVERHEADS
---------------------------------------------------------------------------------------------------------------------------------------------
*/ 
// V4d1  =         ( C1 +   C2               * vol_2      ) / vol_ / C3;                            // [VAR]-dependent updates per loop
// V5d2  =         ( C1 + ( C2 - ttm_ )      * vol_2      ) / vol_ / C3;                            // [VAR]-dependent updates per loop
// Vmd2  =         (           ( ttm_ - C2 ) * vol_2 - C1 ) / vol_ / C3;                            // [VAR]-dependent updates per loop
// 
// V4d1  = U3inv * ( C1 +                 C2 * vol_2      ) / vol_;                                 // fMUL faster than fDIV + a few more fUtilityCONSTs
// V5d2  = U3inv * ( C1 +                 U4 * vol_2      ) / vol_;                                 // fMUL faster than fDIV + a few more fUtilityCONSTs
// Vmd2  = U3inv * (                      U5 * vol_2 - C1 ) / vol_;                                 // fMUL faster than fDIV + a few more fUtilityCONSTs
// 
// ---------------------------------------------------------------------------------------------------
// THIS AVOIDS RE-CALCULATION OF ALL do-while-(CONST)s BY THEIR RE-USE :
//
// if ( option_price < GBlackScholes_WHILEd( vol_,       // <----------[VAR] input_vol,             // GBlackScholes( future_price,
//                                           vol_ * C3,  // <----------[VAR] V1,                    //                strike,
//                                           is_call,    //                  is_call,               //                ttm,
//                                           ttm_,       // do-while-(CONST) ttm_,                  //                risk_rate,
//                                           C1, C2, C3, // do-while-(CONST) C1, C2, C3,            //                b_rate,
//                                           R1, R2      // do-while-(CONST) R1, R2                 //                mid,         // == (high+low)/2,
//                                           )           //                                         //                is_call
//      ) ...                                            //                                         //                )
//
// -------------------------------------------------------------------------------------------------- 
// EVEN BETTER :
//
// if ( option_price < ( is_call ? ( R1 * cdfNormal( 0, 1,  U3inv * ( C1 + C2 * vol_2      ) / vol_ )
//                                 - R2 * cdfNormal( 0, 1,  U3inv * ( C1 + U4 * vol_2      ) / vol_ )
//                                   )
//                               : ( R2 * cdfNormal( 0, 1,  U3inv * (      U5 * vol_2 - C1 ) / vol_ )
//                                 - R1 * cdfNormal( 0, 1, -U3inv * ( C1 + C2 * vol_2      ) / vol_ )
//                                   )
//                       )
// 
//     ) ...
// ________________________________( CALL-OPTIONs )__________________________________________________
if ( option_price < (           ( R1 * cdfNormal( 0, 1,  U3inv * ( C1 + C2 * vol_2      ) / vol_ )
- R2 * cdfNormal( 0, 1,  U3inv * ( C1 + U4 * vol_2      ) / vol_ )
)
)
) {  high = mid;                             // == (high+low)/2;  // LOWER    HI-SIDE BRACKET
} else {   low = mid;                             // == (high+low)/2;  // HEIGHTEN LO-SIDE BRACKET
}
} while ( ( high - low ) > 0.00001 );
return    ( high + low ) / 2; // ________________________________________________________________ JIT/RET
} else {  
do {  mid   = ( high + low ) / 2;       // cheapest  do-while per-loop-[VAR]-update
vol_  = mid + 0.000000000000001; //  cheapest  do-while per-loop-[VAR]-update
vol_2 = vol_ * vol_;            //   cheapest  do-while per-loop-[VAR]-update
// ________________________________( PUT-OPTIONs )___________________________________________________
if ( option_price < (           ( R2 * cdfNormal( 0, 1,  U3inv * (      U5 * vol_2 - C1 ) / vol_ )
- R1 * cdfNormal( 0, 1, -U3inv * ( C1 + C2 * vol_2      ) / vol_ )
)
)
) {  high = mid;                             // == (high+low)/2;  // LOWER    HI-SIDE BRACKET
} else {   low = mid;                             // == (high+low)/2;  // HEIGHTEN LO-SIDE BRACKET
}
} while ( ( high - low ) > 0.00001 );
return    ( high + low ) / 2; // ________________________________________________________________ JIT/RET
}
}

DolphinDB 1.01引入了即时编译(JIT(。JIT 非常适合非矢量化操作,例如隐含波动率的二进制搜索。

若要使用 JIT,只需在用户定义函数之前添加表示法@jit。

最新更新