ブラックショールズ方程式を使用したヨーロピアン・オプションの価格計算

ブラックショールズの単純実装

このコードは、コールおよびプットオプションの価格計算にクローズド形式のソリューションを実装します。

void BlackScholesFormula( int nopt, 
    tfloat r, tfloat sig, const tfloat s0[], const tfloat x[], 
    const tfloat t[], tfloat vcall[], tfloat vput[] ) { 
    tfloat d1, d2; 
    int i; 
    for ( i=0; i<nopt; i++ ) { 
        d1 = ( LOG(s0[i]/x[i]) + (r + HALF*sig*sig)*t[i] ) / ( sig*SQRT(t[i]) ); 
        d2 = ( LOG(s0[i]/x[i]) + (r - HALF*sig*sig)*t[i] ) / ( sig*SQRT(t[i]) ); 

        vcall[i] = s0[i]*CDFNORM(d1) - EXP(-r*t[i])*x[i]*CDFNORM(d2); 
        vput[i] = EXP(-r*t[i])*x[i]*CDFNORM(-d2) - s0[i]*CDFNORM(-d1); 
    } 
}

オプションの数は nopt パラメーターで指定します。tfloat の型は、使用する精度に応じて float または double のいずれかです。同様に、LOG、EXP、SQRT、CDFNORM を数学関数の単精度または倍精度バージョンにマップします。定数 HALF は 0.5f (単精度) または 0.5 (倍精度) のいずれかです。

nopt に加えて、アルゴリズムの入力パラメーターは、s0 (現在の株価)、r (リスク中立レート)、sig (ボラティリティー)、t (満期日)、x (権利行使価格) です。結果は vcall および vput (コールオプションとプットオプションの現在の値) で返されます。

r および sig は価格付けを行うすべてのオプションの定数であり、ほかのパラメーターは浮動小数点値の配列であると仮定しています。vcall および vput パラメーターは出力配列です。

ブラックショールズの最適化実装

void BlackScholesFormula( int nopt, 
    tfloat r, tfloat sig, tfloat s0[], tfloat x[], 
    tfloat t[], tfloat vcall[], tfloat vput[] ) { 
    int i; 
    tfloat a, b, c, y, z, e; 
    tfloat d1, d2, w1, w2; 
    tfloat mr = -r; 
    tfloat sig_sig_two = sig * sig * TWO; 

    for ( i = 0; i < nopt; i++ ) { 
        a = LOG( s0[i] / x[i] ); 
        b = t[i] * mr; 
        z = t[i] * sig_sig_two; 
        c = QUARTER * z; 
        e = EXP ( b ); 
        y = INVSQRT( z ); 

        w1 = ( a - b + c ) * y; 
        w2 = ( a - b - c ) * y; 
        d1 = ERF( w1 ); 
        d2 = ERF( w2 ); 
        d1 = HALF + HALF*d1; 
        d2 = HALF + HALF*d2; 

        vcall[i] = s0[i]*d1 - x[i]*e*d2; 
        vput[i] = vcall[i] - s0[i] + x[i]*e; 
    } 
}

このコードで、INVSQRT(x) は精度に応じて 1.0/sqrt(x) または 1.0f/sqrtf(x) のいずれかです。TWO および QUARTER はそれぞれ、浮動小数点定数 2 と 0.25 です。

ブラックショールズのイン oneMKL VM 実装

void BlackScholesFormula_MKL( int nopt, 
    tfloat r, tfloat sig, tfloat * s0, tfloat * x, 
    tfloat * t, tfloat * vcall, tfloat * vput ) { 
    int i; 
    tfloat mr = -r; 
    tfloat sig_sig_two = sig * sig * TWO; 

    #pragma omp parallel for \ 
        shared(s0, x, t, vcall, vput, mr, sig_sig_two, nopt) \ 
        default(none) 
    for ( i = 0; i < nopt; i+= NBUF ) { 
        int j; 
        tfloat *a, *b, *c, *y, *z, *e; 
        tfloat *d1, *d2, *w1, *w2; 
        __declspec(align(ALIGN_FACTOR)) tfloat Buffer[NBUF*4]; 
        // nopt がNBUF の倍数でなかった場合 
        // ループの最後の反復のベクトル長を計算 
        #define MY_MIN(x, y) ((x) < (y)) ? (x) : (y) 
        int nbuf = MY_MIN(NBUF, nopt - i); 

        a = Buffer + NBUF*0; w1 = a; d1 = w1; 
        c = Buffer + NBUF*1; w2 = c; d2 = w2; 
        b = Buffer + NBUF*2; e = b; 
        z = Buffer + NBUF*3; y = z; 

        // 各スレッドにVML の精度を設定 
        vmlSetMode( VML_ACC ); 

        VDIV(nbuf, s0 + i, x + i, a); 
        VLOG(nbuf, a, a); 

        #pragma simd 
        for ( j = 0; j < nbuf; j++ ) { 
            b[j] = t[i + j] * mr; 
            a[j] = a[j] - b[j]; 
            z[j] = t[i + j] * sig_sig_two; 
            c[j] = QUARTER * z[j]; 
        } 

        VINVSQRT(nbuf, z, y); 
        VEXP(nbuf, b, e); 

        #pragma simd 
        for ( j = 0; j < nbuf; j++ ) { 
            tfloat aj = a[j]; 
            tfloat cj = c[j]; 
            w1[j] = ( aj + cj ) * y[j]; 
            w2[j] = ( aj - cj ) * y[j]; 
        } 

        VERF(nbuf, w1, d1); 
        VERF(nbuf, w2, d2); 

        #pragma simd 
        for ( j = 0; j < nbuf; j++ ) { 
            d1[j] = HALF + HALF*d1[j]; 
            d2[j] = HALF + HALF*d2[j]; 
            vcall[i+j] = s0[i+j]*d1[j] - x[i+j]*e[j]*d2[j]; 
            vput[i+j] = vcall[i+j] - s0[i+j] + x[i+j]*e[j]; 
        } 
    } 
}

同等の精度では、計算により制限される問題 (データが L2 キャッシュに収まる) の場合、oneMKL VM はインテル® コンパイラーと SVML ベースのソリューションよりも 30-50% 優れたパフォーマンスを実現します。この場合、キャッシュ読み取り/書き込み操作のレイテンシーは計算によりマスクされます。問題サイズが大きくなりメモリー帯域幅により制限されるようになると、メモリー使用量の最適化はさらに重要となり、中間バッファーを利用するインテルの VM ソリューションは SVML を利用したバッファリングなしの 1 パス・ソリューションに対する優位性を失います。