Parallel Computing Toolbox

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

この例では、金融オプションの価格をモンテカルロ法を使用して GPU で計算する方法を示します。エキゾチック オプションの 3 つの単純なタイプを例として使用しますが、同様の方法でもっと複雑なオプションの価格を設定できます。

この例は関数であるため、補助関数をその内部に入れ子にすることができます。

function paralleldemo_gpu_optionpricing

この例では実行に長時間かかるカーネルを使用するため、GPU でのカーネル実行がタイムアウト可能な場合は実行できません。タイムアウトが有効なのは、通常は、選択された GPU がディスプレイの駆動も行っている場合に限られます。

dev = gpuDevice();
if dev.KernelExecutionTimeout
    error( 'pctexample:gpuoptionpricing:KernelTimeout', ...
        ['This example cannot run if kernel execution on the GPU can ', ...
        'time-out.'] );
end

株価の変化

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

$${/rm d}S = S /times
  /left[
     (r-d){/rm d}t + /sigma /epsilon /sqrt{{/rm d}t}
  /right]$$

ここで、 $S$ は株価、 $r$ は無リスク金利、 $d$ は株の年次配当利回り、 $/sigma$ は価格の予想変動率、 $/epsilon$ はガウスのホワイト ノイズ過程を表します。 $(S+/Delta S)/S$ が対数正規分布していると仮定すると、これは次のように離散化できます。

$$S_{t+1} = S_{t} /times /exp{
  /left[
    /left(r - d - /frac{1}{2}/sigma^2/right)/Delta t
    + /sigma /epsilon /sqrt{ /Delta t }
  /right]
} $$

例として、毎年 1% の配当を生む $100 の株を使用します。政府の基準貸し付け利率が 0.5% であるとします。ほぼ毎日サンプリングした 2 年間の時間ウィンドウを調べます。市場の変動は、毎年 20% であると仮定します。

stockPrice   = 100;   % Stock price starts at $100.
dividend     = 0.01;  % 1% annual dividend yield.
riskFreeRate = 0.005; % 0.5 percent.
timeToExpiry = 2;     % Lifetime of the option in years.
sampleRate   = 1/250; % Assume 250 working days per year.
volatility   = 0.20;  % 20% volatility.

乱数発生器をリセットして、結果が必ず繰り返されるようにします。

seed = 1234;
rng( seed );              % Reset the CPU random number generator.
parallel.gpu.rng( seed ); % Reset the GPU random number generator.

これで、長期にわたるループにより株価の変動経路をシミュレートできます。

price = stockPrice;
time = 0;
hold on;
while time < timeToExpiry
    time = time + sampleRate;
    drift = (riskFreeRate - dividend - volatility*volatility/2)*sampleRate;
    perturbation = volatility*sqrt( sampleRate )*randn();
    price = price*exp(drift + perturbation);
    plot( time, price, '.' );
end
axis tight;
grid on;
xlabel( 'Time (years)' );
ylabel( 'Stock price ($)' );

GPU での実行

株価シミュレーションを GPU で実行するには、まずシミュレーションのループを補助関数内に配置する必要があります。

    function finalStockPrice = simulateStockPrice(S,r,d,v,T,dT)
        t = 0;
        while t < T
            t = t + dT;
            dr = (r - d - v*v/2)*dT;
            pert = v*sqrt( dT )*randn();
            S = S*exp(dr + pert);
        end
        finalStockPrice = S;
    end

こうすると、arrayfun を使用してこれを何千回も呼び出すことができます。計算が GPU で行われるように、入力価格は、シミュレーションごとに 1 つの要素をもつ GPU ベクトルにします。GPU での計算時間を正確に測定するため、関数 gputimeit を使用します。

% Create the input data.
N = 1000000;
startStockPrices = stockPrice*gpuArray.ones(N,1);

% Run the simulations.
finalStockPrices = arrayfun( @simulateStockPrice, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    timeToExpiry, sampleRate );
meanFinalPrice = mean(finalStockPrices);

% Measure the execution time of the function on the GPU using gputimeit.
% This requires us to store the |arrayfun| call in a function handle.
functionToTime = @() arrayfun(@simulateStockPrice, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    timeToExpiry, sampleRate );
timeTaken = gputimeit(functionToTime);

fprintf( 'Calculated average price of $%1.4f in %1.3f secs.\n', ...
    meanFinalPrice, timeTaken );

clf;
hist( finalStockPrices, 100 );
xlabel( 'Stock price ($)' )
ylabel( 'Frequency' )
grid on;
Calculated average price of $98.9563 in 0.275 secs.

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

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

    function optionPrice = asianCallOption(S,r,d,v,x,T,dT)
        t = 0;
        cumulativePrice = 0;
        while t < T
            t = t + dT;
            dr = (r - d - v*v/2)*dT;
            pert = v*sqrt( dT )*randn();
            S = S*exp(dr + pert);
            cumulativePrice = cumulativePrice + S;
        end
        numSteps = (T/dT);
        meanPrice = cumulativePrice / numSteps;
        % Express the final price in today's money.
        optionPrice = exp(-r*T) * max(0, meanPrice - x);
    end

ここでも GPU を使用して、arrayfun により何千ものシミュレーション パスを実行します。各シミュレーション パスはオプション価格の独立した見積もりを与えるため、平均値を結果と見なします。

strike = 95;   % Strike price for the option ($).

optionPrices = arrayfun( @asianCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, strike, ...
    timeToExpiry, sampleRate );
meanOptionPrice = mean(optionPrices);

% Measure the execution time on the GPU and show the results.
functionToTime = @() arrayfun( @asianCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, strike, ...
    timeToExpiry, sampleRate );
timeTaken = gputimeit(functionToTime);

fprintf( 'Calculated average price of $%1.4f in %1.3f secs.\n', ...
    meanOptionPrice, timeTaken );
Calculated average price of $8.7210 in 0.279 secs.

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

この例ではヨーロピアンスタイルのルックバック オプションを使用します。その支払いは、オプションの有効期間における最低株価と最終株価の差額になります。

    function optionPrice = euroLookbackCallOption(S,r,d,v,T,dT)
        t = 0;
        minPrice = S;
        while t < T
            t = t + dT;
            dr = (r - d - v*v/2)*dT;
            pert = v*sqrt( dT )*randn();
            S = S*exp(dr + pert);
            if S<minPrice
                minPrice = S;
            end
        end
        % Express the final price in today's money.
        optionPrice = exp(-r*T) * max(0, S - minPrice);
    end

この場合、オプションの行使価格は最低株価であることに注意してください。最終株価は必ず最低株価以上となるため、オプションは常に行使され、本来の意味で「オプション」とは言えません。

optionPrices = arrayfun( @euroLookbackCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    timeToExpiry, sampleRate );
meanOptionPrice = mean(optionPrices);

% Measure the execution time on the GPU and show the results.
functionToTime = @() arrayfun( @euroLookbackCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    timeToExpiry, sampleRate );
timeTaken = gputimeit(functionToTime);

fprintf( 'Calculated average price of $%1.4f in %1.3f secs.\n', ...
    meanOptionPrice, timeTaken );
Calculated average price of $19.2711 in 0.277 secs.

バリア オプションの価格

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

    function optionPrice = upAndOutCallOption(S,r,d,v,x,b,T,dT)
        t = 0;
        while (t < T) && (S < b)
            t = t + dT;
            dr = (r - d - v*v/2)*dT;
            pert = v*sqrt( dT )*randn();
            S = S*exp(dr + pert);
        end
        if S<b
            % Within barrier, so price as for a European option.
            optionPrice = exp(-r*T) * max(0, S - x);
        else
            % Hit the barrier, so the option is withdrawn.
            optionPrice = 0;
        end
    end

ここでは、オプションの行使価格とそれが無効となるバリア価格の両方を提供しなければならないことに注意してください。

strike  = 95;   % Strike price for the option ($).
barrier = 150;  % Barrier price for the option ($).

optionPrices = arrayfun( @upAndOutCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    strike, barrier, ...
    timeToExpiry, sampleRate );
meanOptionPrice = mean(optionPrices);

% Measure the execution time on the GPU and show the results.
functionToTime = @() arrayfun( @upAndOutCallOption, ...
    startStockPrices, riskFreeRate, dividend, volatility, ...
    strike, barrier, ...
    timeToExpiry, sampleRate );
timeTaken = gputimeit(functionToTime);

fprintf( 'Calculated average price of $%1.4f in %1.3f secs.\n', ...
    meanOptionPrice, timeTaken );
Calculated average price of $6.8166 in 0.279 secs.
end