%%
% 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 mc2c

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('CPU (arrayfun): value = %f, execution time = %f seconds \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('GPU (arrayfun): value = %f, execution time = %f seconds \n',ave,toc);

%
% now do it in a vectorised form
%

tic
vals = path2d_v(S1,S2,r,sig1,sig2,rho,T,N,K, M);
ave  = mean(vals);

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

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

wait(gpu);
fprintf('GPU (vector 1): value = %f, execution time = %f seconds \n',ave,toc);

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

wait(gpu);
fprintf('GPU (vector 2): value = %f, execution time = %f seconds \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

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

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

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

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

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

end


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

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

for n = 1:N
  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + 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


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

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

%
% loop unrolling reduces CUDA kernel launch overhead?
%

for n1 = 1:8:N
  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + rho*dW1;  % correlated Brownian increments
  S1  = S1.*(1 + r*dt + sig1*dW1);
  S2  = S2.*(1 + r*dt + sig2*dW2);

  dW1 = con1*gpuArray.randn(1,M);            % turn uncorrelated Normals into
  dW2 = con2*gpuArray.randn(1,M) + 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
