Main Content

モンテカルロ シミュレーションでの GPU arrayfun の使用

この例では、金融オプションの価格をモンテカルロ法を使用して GPU で計算する方法を示します。

この例ではエキゾチック オプションの 3 つの単純なタイプを使用しますが、同様の方法でもっと複雑なオプションの価格を設定できます。この例では、モンテカルロ シミュレーションを CPU で実行する場合と、GPU でarrayfunを使用して実行する場合の実行時間を比較します。

株価の変化

ここでは、無リスク金利、配当利回り (該当する場合)、市場の予想変動率に関連した対数正規分布に従って価格が変化するものと仮定します。さらに、こうした量はすべて、オプションの有効期間中は固定されていると仮定します。これらの仮定によって、価格に関する次のような確率微分方程式が導かれます。

dS=S×[(r-d)dt+σϵdt],

ここで、S は株価、r は無リスク金利、d は株の年間配当利回り、σ は価格のボラティリティ、ϵ はガウスのホワイト ノイズ過程です。(S+ΔS)/S が対数正規分布であると仮定すると、この微分方程式は離散化でき、次の方程式が得られます。

St+1=St×exp[(r-d-12σ2)Δt+σϵΔt].

以上の仮定の下で、$100 の株を使用して 2 年間の時間枠を調べます。

  • この株は毎年 1% の配当を生みます。

  • 政府の無リスク金利は 0.5% です。

  • 価格は毎日サンプリングされ、稼働日は年間 250 日です。

  • 市場の変動は毎年 20% です。

stockPrice = 100;
timeToExpiry = 2;
dividend = 0.01;
riskFreeRate = 0.005;
sampleRate = 1/250;
volatility = 0.20;

必ず予測可能な結果が得られるように、CPU と GPU の乱数発生器のシードを設定します。

seed = 1234;
rng(seed);
gpurng(seed);

長期にわたる株価の変動経路をシミュレートし、結果をプロットします。

price = stockPrice;
time = 0;
h = animatedline(Marker=".");

while time < timeToExpiry
    time = time + sampleRate;
    drift = (riskFreeRate - dividend - volatility*volatility/2)*sampleRate;
    perturbation = volatility*sqrt(sampleRate)*randn;
    price = price*exp(drift + perturbation);
    addpoints(h,time,price);
end

grid on
axis tight
xlabel("Time (years)")
ylabel("Stock Price ($)")

CPU と GPU での実行時間の測定

この例の最後に示している関数 simulateStockPrice は、前のセクションで説明した離散化された微分方程式を使用して株価をシミュレートします。

株価の 100,000 個のモンテカルロ シミュレーションを実行するための入力データを準備します。

N = 100000;
startStockPrices = stockPrice*ones(N,1);

CPU における 100,000 個のシミュレーションの実行時間を測定します。

tic
finalStockPricesCPU = zeros(N,1);
for i = 1:N
    finalStockPricesCPU(i) = simulateStockPrice(startStockPrices(i), ...
        riskFreeRate,dividend,volatility, ...
        timeToExpiry,sampleRate);
end
timeCPU = toc;

各シミュレーションはオプション価格の独立した見積もりを与えるため、平均値を結果と見なします。

fprintf("Calculated average price of $%1.4f on the CPU in %1.3f secs.\n", ...
    mean(finalStockPricesCPU),timeCPU);
Calculated average price of $99.0857 on the CPU in 2.206 secs.

シミュレーションを GPU で実行するには、gpuArrayオブジェクトを作成して入力データを GPU 上で準備します。

gpuStartStockPrices = gpuArray(startStockPrices);

GPU 配列と関数ハンドルを入力として arrayfun を呼び出すと、arrayfun は、指定した関数を配列の各要素に適用します。この動作は、各開始株価のループ処理が不要であることを意味します。GPU の関数 arrayfun を使って要素単位の MATLAB® 関数をカスタム CUDA® カーネルに変換すると、処理実行のオーバーヘッドを削減できます。

arrayfun を使用して関数 simulateStockPrice を実行し、gputimeitを使用して GPU における 100,000 個のシミュレーションの実行時間を測定します。

finalStockPricesGPU = arrayfun(@simulateStockPrice, ...
        gpuStartStockPrices,riskFreeRate,dividend,volatility, ...
        timeToExpiry,sampleRate);

timeGPU = gputimeit(@() arrayfun(@simulateStockPrice, ...
        gpuStartStockPrices,riskFreeRate,dividend,volatility, ...
        timeToExpiry,sampleRate));

fprintf("Calculated average price of $%1.4f on the GPU in %1.3f secs.\n", ...
    mean(finalStockPricesGPU),timeGPU);
Calculated average price of $99.0442 on the GPU in 0.023 secs.

GPU におけるモンテカルロ シミュレーションの結果をヒストグラムにプロットします。

histogram(finalStockPricesGPU,100);
xlabel("Stock Price ($)")
ylabel("Frequency")
grid on

アジアン オプションの価格

オプションの有効期間における株価の算術平均を基に、ヨーロピアン タイプのアジアン オプションを使用します。関数 asianCallOption は、シミュレーション中に価格を累積することで、平均価格を計算します。コール オプションの場合、この関数は、平均価格が行使価格を超えている場合にオプションを行使します。支払いは平均価格と行使価格の間の差額となります。この例の最後に示している asianCallOption を使用して、アジアン コール オプションをシミュレートします。

オプションの行使価格を $95 に設定します。

strike = 95;

arrayfun を使用して CPU と GPU における 100,000 個のシミュレーションの実行時間を測定し、結果を表示します。

tic
optionPricesCPU = zeros(N,1);
for i=1:N
    optionPricesCPU(i) = asianCallOption(startStockPrices(i), ...
        riskFreeRate,dividend,volatility,strike, ...
        timeToExpiry,sampleRate);
end
timeAsianOptionCPU = toc;

fprintf("Calculated average price of $%1.4f on the CPU in %1.3f secs.\n", ...
    mean(optionPricesCPU),timeAsianOptionCPU);
Calculated average price of $8.6733 on the CPU in 2.146 secs.
optionPricesGPU = arrayfun( @asianCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility,strike, ...
    timeToExpiry,sampleRate );

timeAsianOptionGPU = gputimeit(@() arrayfun( @asianCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility,strike, ...
    timeToExpiry,sampleRate));

fprintf("Calculated average price of $%1.4f on the GPU in %1.3f secs.\n", ...
    mean(optionPricesGPU),timeAsianOptionGPU );
Calculated average price of $8.7448 on the GPU in 0.023 secs.

ルックバック オプションの価格

ヨーロピアンスタイルのルックバック オプションを使用します。その支払いは、オプションの有効期間における最低株価と最終株価の差額になります。オプションの行使価格は最低株価です。最終株価は必ず最低株価以上となるため、オプションは常に行使され、本来の意味でオプションとは言えません。この例の最後に示している lookbackCallOption を使用して、ヨーロピアンスタイルのルックバック コール オプションをシミュレートします。

arrayfun を使用して CPU と GPU の両方における 100,000 個のシミュレーションの実行時間を測定し、結果を表示します。

tic
optionPricesCPU = zeros(N,1);
for i=1:N
    optionPricesCPU(i) = lookbackCallOption(startStockPrices(i), ...
        riskFreeRate,dividend,volatility, ...
        timeToExpiry,sampleRate);
end
timeLookbackOptionCPU = toc;

fprintf("Calculated average price of $%1.4f on the CPU in %1.3f secs.\n", ...
    mean(optionPricesCPU),timeLookbackOptionCPU);
Calculated average price of $19.2456 on the CPU in 2.201 secs.
optionPricesGPU = arrayfun(@lookbackCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility, ...
    timeToExpiry,sampleRate);

timeLookbackOptionGPU = gputimeit(@() arrayfun(@lookbackCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility, ...
    timeToExpiry,sampleRate));

fprintf("Calculated average price of $%1.4f on the GPU in %1.3f secs.\n", ...
    mean(optionPricesGPU),timeLookbackOptionGPU);
Calculated average price of $19.3893 on the GPU in 0.021 secs.

バリア オプションの価格

株価がバリア レベルに到達すると無効になるアップ アンド アウト バリア オプションを使用します。株価がバリア レベル未満に留まる場合は、最終株価を通常のヨーロピアン コール オプションの計算に使用します。この例の最後に示している関数 upAndOutCallOption を使用して、アップ アンド アウト バリア コール オプションをシミュレートします。

オプションの行使価格とそれが無効となるバリア価格を設定します。行使価格を $95、バリア価格を $150 にします。

strike  = 95;
barrier = 150;

arrayfun を使用して CPU と GPU における 100,000 個のシミュレーションの実行時間を測定し、結果を表示します。

tic
optionPricesCPU = zeros(N,1);
for i=1:N
    optionPricesCPU(i) = upAndOutCallOption(startStockPrices(i), ...
        riskFreeRate,dividend,volatility,strike, ...
        barrier,timeToExpiry,sampleRate);
end
timeBarrierOptionCPU = toc;

fprintf("Calculated average price of $%1.4f on the CPU in %1.3f secs.\n", ...
    mean(optionPricesCPU),timeBarrierOptionCPU);
Calculated average price of $6.8327 on the CPU in 2.074 secs.
optionPricesGPU = arrayfun(@upAndOutCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility,strike, ...
    barrier,timeToExpiry,sampleRate);

timeBarrierOptionGPU = gputimeit(@() arrayfun(@upAndOutCallOption, ...
    gpuStartStockPrices,riskFreeRate,dividend,volatility,strike, ...
    barrier,timeToExpiry,sampleRate));

fprintf("Calculated average price of $%1.4f on the GPU in %1.3f secs.\n", ...
    mean(optionPricesGPU),timeBarrierOptionGPU);
Calculated average price of $6.7834 on the GPU in 0.021 secs.

結果の比較

シミュレーションごとに、GPU の実行時間に対する CPU の実行時間の比率を計算します。

ratio = [timeCPU/timeGPU timeAsianOptionCPU/timeAsianOptionGPU ...
    timeLookbackOptionCPU/timeLookbackOptionGPU timeBarrierOptionCPU/timeBarrierOptionGPU]
ratio = 1×4

   94.2557   94.6009  104.1725   97.5490

結果を可視化するために、各シミュレーションの実行時間の比率をプロットします。

bar(categorical(["Stock Price" "Asian Call Option" "Lookback Option" "Barrier Option"]), ...
    ratio)
ylabel("Ratio of CPU to GPU Execution Time")

この例では、arrayfun を使用してシミュレーションを GPU で実行した場合の方が、シミュレーションを CPU で実行した場合よりも大幅に高速になっています。

この例で説明した手法を独自のコードに適用した場合のパフォーマンスの改善は、ハードウェアや実行するコードに大きく依存します。

サポート関数

株価の変化のシミュレーション関数

関数 simulateStockPrice は、モンテカルロ シミュレーションを実行して最終株価を特定します。この計算では、無リスク金利、配当利回り (該当する場合)、市場の予想変動率に関連した対数正規分布に従って価格が変化するものと仮定します。

関数 simulateStockPrice は、初期株価、無リスク金利、配当率、市場の変動、合計時間枠、およびサンプル レートを入力として受け取ります。

function finalStockPrice = simulateStockPrice(price,riskFreeRate,dividend,volatility,T,dT)
t = 0;
while t < T
    t = t + dT;
    drift = (riskFreeRate - dividend - volatility*volatility/2)*dT;
    perturbation = volatility*sqrt(dT)*randn;
    price = price.*exp(drift + perturbation);
end

finalStockPrice = price;
end

アジアン コール オプション関数

関数 asianCallOption は、モンテカルロ シミュレーションを実行してアジアン コール オプションの価格を特定します。この計算では、無リスク金利、配当利回り (該当する場合)、市場の予想変動率に関連した対数正規分布に従って価格が変化するものと仮定します。この関数は、シミュレーション中に価格を累積することで、平均価格を計算します。コール オプションの場合、この関数は、平均価格が行使価格を超えている場合にオプションを行使します。支払いは平均価格と行使価格の間の差額となります。

関数 asianCallOption は、初期株価、無リスク金利、配当率、市場の変動、行使価格、合計時間枠、およびサンプル レートを入力として受け取ります。

function optionPrice = asianCallOption(price,riskFreeRate,dividend,volatility,strike,T,dT)
t = 0;
cumulativePrice = 0;
while t < T
    t = t + dT;
    dr = (riskFreeRate - dividend - volatility*volatility/2)*dT;
    pert = volatility*sqrt(dT)*randn;
    price = price*exp(dr + pert);
    cumulativePrice = cumulativePrice + price;
end
numSteps = (T/dT);
meanPrice = cumulativePrice/numSteps;

% Express the final price in today's money
optionPrice = exp(-riskFreeRate*T)*max(0,meanPrice - strike);
end

ルックバック オプション関数

関数 lookbackCallOption は、モンテカルロ シミュレーションを実行してヨーロピアンスタイルのルックバック オプションを特定します。その支払いは、オプションの有効期間における最低株価と最終株価の差額になります。オプションの行使価格は最低株価です。最終株価は必ず最低株価以上となるため、オプションは常に行使され、本来の意味で「オプション」とは言えません。

関数 lookbackCallOption は、初期株価、無リスク金利、配当率、市場の変動、合計時間枠、およびサンプル レートを入力として受け取ります。

function optionPrice = lookbackCallOption(price,riskFreeRate,dividend,volatility,T,dT)
t = 0;
minPrice = price;
while t < T
    t = t + dT;
    dr = (riskFreeRate - dividend - volatility*volatility/2)*dT;
    pert = volatility*sqrt(dT)*randn;
    price = price*exp(dr + pert);
    if price < minPrice
        minPrice = price;
    end
end

% Express the final price in today's money
optionPrice = exp(-riskFreeRate*T)*max(0,price - minPrice);
end

バリア オプション関数

関数 upAndOutCallOption は、モンテカルロ シミュレーションを実行してアップ アンド アウト バリア コール オプションの価格を特定します。株価がバリア レベル未満に留まる場合、この関数は最終株価を通常のヨーロピアン コール オプションの計算に使用します。

関数 upAndOutCallOption は、初期株価、無リスク金利、配当率、市場の変動、行使価格、バリア価格、合計時間枠、およびサンプル レートを入力として受け取ります。

function optionPrice = upAndOutCallOption(price,riskFreeRate,dividend,volatility,strike,barrier,T,dT)
t = 0;
while (t < T) && (price < barrier)
    t = t + dT;
    dr = (riskFreeRate - dividend - volatility*volatility/2)*dT;
    pert = volatility*sqrt(dT)*randn;
    price = price*exp(dr + pert);
end

if price<barrier
    % Within barrier, so price as for a European option
    optionPrice = exp(-riskFreeRate*T)*max(0,price - strike);
else
    % Hit the barrier, so the option is withdrawn
    optionPrice = 0;
end
end

参考

| |

関連するトピック