Использование графического процессора ARRAYFUN для моделирования Монте-Карло

Открыть сценарий

В этом примере показано, как цены для финансовых опционов могут быть рассчитаны на GPU с использованием методов Монте-Карло. В качестве примеров используются три простых типа экзотических опций, но более сложные опции могут быть оценены аналогичным образом.

Этот пример является функцией, позволяющей вложить в нее вспомогательные элементы.

function paralleldemo_gpu_optionpricing

В этом примере используются давно работающие ядра, поэтому их выполнение невозможно, если выполнение ядра на 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]
}$

В качестве примера давайте используем $100 акций, которые приносят 1% дивидендов каждый год. Предполагается, что процентная ставка центрального правительства составит 0,5%. Мы исследуем двухлетнее временное окно, отобранное примерно ежедневно. Волатильность рынка предполагается на уровне 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.
gpurng( 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 сначала необходимо поместить цикл моделирования в вспомогательную функцию:

    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 мы делаем входные цены вектором GPU с одним элементом на моделирование. Для точного измерения времени расчета на GPU мы используем gputimeit функция.

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

% 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.283 secs.

Ценообразование азиатского опциона

В качестве примера давайте используем европейский азиатский опцион на основе среднего арифметического цены акций в течение срока действия опциона. Среднюю цену можно рассчитать путем накопления цены во время моделирования. Для колл-опциона опцион реализуется, если средняя цена выше страйка, а выплата представляет собой разницу между двумя:

    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

Снова мы используем графический процессор для запуска тысяч путей моделирования с использованием 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.287 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.286 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.289 secs.

end

Документация