Используя графический процессор 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
Эволюция курса акций
Мы принимаем, что цены развиваются согласно логарифмически нормальному распределению, связанному с безрисковой процентной ставкой, дивидендной доходностью (если таковые имеются), и энергозависимость на рынке. Все эти количества приняты зафиксированные за время жизни опции. Это дает следующее стохастическое дифференциальное уравнение за цену:
![$${\rm d}S = S \times
 \left[
 (r-d){\rm d}t + \sigma \epsilon \sqrt{{\rm d}t}
 \right]$$](../../examples/parallel/win64/paralleldemo_gpu_optionpricing_eq17284130511455291481.png)
где
курс акций,
является безрисковой процентной ставкой,
является ежегодной дивидендной доходностью запаса,
является энергозависимостью цены и
представляет Гауссов бело-шумовой процесс. Принятие этого
логарифмически нормально распределяется, это может быть дискретизировано к:
![$$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]
} $$](../../examples/parallel/win64/paralleldemo_gpu_optionpricing_eq14584843645873655464.png)
Как пример позволяют нам использовать 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