Использование GPU ARRAYFUN для симуляций Монте-Карло
В этом примере показано, как цены для финансовых опций могут быть рассчитаны на графическом процессоре с помощью методов Monte-Carlo. В качестве примеров используются три простых типа экзотической опции, но аналогичным образом можно оценить и более сложные опции.
Этот пример является функцией, так что помощники могут быть вложены в него.
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%. Мы исследуем двухлетний временной интервал, отбираемый примерно ежедневно. Волатильность рынка принята 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 ($)' );
Выполнение на графическом процессоре
Чтобы запустить симуляции цен на акции на графическом процессоре, нам сначала нужно поместить цикл симуляции в функцию helper:
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 с одним элементом за симуляцию. Чтобы точно измерить время вычисления на графическом процессоре, мы используем 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.
Ценообразование азиатской опции
В качестве примера давайте будем использовать европейский азиатский вариант, основанный на среднеарифметическом значении цены акции в течение срока службы опции. Мы можем вычислить среднюю цену, накопив цену во время симуляции. Для опции call, опция реализуется, если средняя цена выше страйка, и выплата является различием между ними:
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