%%
% This is my example of Monte Carlo pricing based on 2D path simulations.
% It is based on the MATLAB demo application paralleldemo_gpu_optionpricing
%
% Some notes on the application are available here:
% http://people.maths.ox.ac.uk/gilesm/cuda/prac2/MC_notes.pdf
%

function mc2

S1_0 = 100;
S2_0 = 100;
r    = 0.05;
sig1 = 0.2;
sig2 = 0.3;
rho  = 0.2;
T    = 1;
K    = 100;

% reset the random number generators to ensure repeatable results

gpu = gpuDevice;

seed = 1234;
rng(seed); parallel.gpu.rng(seed);

M = 1e6;  % number of paths
N = 128;  % number of timesteps

%
% first do it on the CPU -- each input element has to be a vector
% it appears the execution is not multi-threaded
%

tic
o = ones(1,M);
vals = arrayfun(@path2d, S1_0*o,S2_0*o,r*o,sig1*o,sig2*o,rho*o,T*o,N*o,K*o);
ave  = mean(vals);

fprintf('ave value = %f, execution time = %f \n',ave,toc);

%
% then do it on the GPU -- scalar input elements automatically replicated
%

tic
S1   = S1_0*gpuArray.ones(1,M);
S2   = S2_0*gpuArray.ones(1,M);
vals = arrayfun(@path2d, S1,S2,r,sig1,sig2,rho,T,N,K);
ave  = mean(vals);
wait(gpu);

fprintf('ave value = %f, execution time = %f \n',ave,toc);

end


function val = path2d(S1,S2,r,sig1,sig2,rho,T,N,K)

dt   = T/N;
con1 = sqrt(dt);
con2 = sqrt((1-rho^2)*dt);

for n = 1:N
  dW1 = con1*randn();            % turn uncorrelated Normals into
  dW2 = con2*randn() + rho*dW1;  % correlated Brownian increments

  S1  = S1*(1 + r*dt + sig1*dW1);
  S2  = S2*(1 + r*dt + sig2*dW2);
end

val = exp(-r*T) * max( min(S1,S2)-K, 0 );   % call option value

end
