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

Открыть скрипт

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

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

function paralleldemo_gpu_optionpricing

Продолжительные ядра использования этого примера, так не может запуститься, если выполнение ядра на графическом процессоре может тайм-аут. Тайм-аут обычно только активен, если выбранный графический процессор также управляет отображением.

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

Эволюция курса акций

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

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

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

Работа графического процессора

Чтобы запустить симуляции курса акций на графическом процессоре, мы сначала должны поместить цикл симуляции в функции помощника:

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

Оценка опции Lookback

Для этого примера мы используем европейский стиль lookback опция, выплата которой является различием между минимальным курсом акций и итоговым курсом акций за время жизни опции.

    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
Документация Parallel Computing Toolbox
Поддержка
Для просмотра документации необходимо авторизоваться на сайте