【CUDA】Monte Carlo Simulationを用いたOption Pricing

別記事で「Cooperative Group」を用いたReductionについて記事を書いたけど、肝心のMonteCarloについては触れていなかったので簡単にメモ。
gomiwo.hatenablog.com

モンテカルロでのプライシング理論はググれば沢山情報が出てくるのでそちらを参照。
簡単に書くと、UndelyingのSDE使って満期までのパスを生成→パスごとにペイオフを算出→パス本数で割る、という流れ。

以下はCUDAのサンプルに入っているMonteCarlo Simulationのサンプルコードを抜粋。

・求めたいOprion Priceが複数ある状況で、１つのOptionにつき１Block割り当て、１つの株価パスにつき1Thread割り当てている。
・「cooperative_group.h」を使ってBlock内のThreadを分割。Reductionで使用。
・Call Option Priceの配列はShared Memoryに格納
・株価パスの生成は#pragma unroll 8でループ展開することで速度向上

static __global__ void MonteCarloOneBlockPerOption(
    curandState * __restrict rngStates,
    const __TOptionData * __restrict d_OptionData,
    __TOptionValue * __restrict d_CallValue,
    int pathN,
    int optionN)
{
    // Handle to thread block group
    cg::thread_block cta = cg::this_thread_block();
    cg::thread_block_tile<32> tile32 = cg::tiled_partition<32>(cta); //Thread Groupを作成

    const int SUM_N = THREAD_N;
    __shared__ real s_SumCall[SUM_N];
    __shared__ real s_Sum2Call[SUM_N];

    // determine global thread id
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    // Copy random number state to local memory for efficiency
    curandState localState = rngStates[tid];
    for(int optionIndex = blockIdx.x; optionIndex < optionN; optionIndex += gridDim.x) //Block1つにつきOption Price1つを計算
    {
        const real        S = d_OptionData[optionIndex].S;
        const real        X = d_OptionData[optionIndex].X;
        const real    MuByT = d_OptionData[optionIndex].MuByT;
        const real VBySqrtT = d_OptionData[optionIndex].VBySqrtT;

        //Cycle through the entire samples array:
        //derive end stock price for each path
        //accumulate partial integrals into intermediate shared memory buffer
        for (int iSum = threadIdx.x; iSum < SUM_N; iSum += blockDim.x) //Thread１つにつきStockのPath1本を生成
        {
            __TOptionValue sumCall = {0, 0};

            #pragma unroll 8 //ループ展開で速度向上
            for (int i = iSum; i < pathN; i += SUM_N)
            {
                real              r = curand_normal(&localState); //正規乱数の生成
                real      callValue = endCallValue(S, X, r, MuByT, VBySqrtT); //PayOffの計算
                sumCall.Expected   += callValue;
                sumCall.Confidence += callValue * callValue;
            }

            s_SumCall[iSum]  = sumCall.Expected;
            s_Sum2Call[iSum] = sumCall.Confidence;
        }

        //Reduce shared memory accumulators
        //and write final result to global memory
        cg::sync(cta);
        sumReduce<real, SUM_N, THREAD_N>(s_SumCall, s_Sum2Call, cta, tile32, &d_CallValue[optionIndex]); //Reductionの実行
    }
}

__device__ inline float endCallValue(float S, float X, float r, float MuByT, float VBySqrtT)
{
    float callValue = S * __expf(MuByT + VBySqrtT * r) - X;
    return (callValue > 0.0F) ? callValue : 0.0F;
}