%GPU_Jacket_corr1 - A first attempt at using Jacket.
%This is a demonstration of how not to do it. See GPU_Jacket_corr3 for a
%MUCH better version
%% Correlated asset information
CurrentPrice = gdouble([78 102]);       %Initial Prices of the two stocks
Corr = [1 0.4; 0.4 1];                  %Correlation Matrix
T = 500;                                %Number of days to simulate = 2years = 500days 
Div=gdouble([0.01 0.01]);               %Dividend
Vol=gdouble([0.2 0.3]);                 %Volatility

%% Market Information
r = gdouble(0.03);                      %Risk-free rate

%% Simulation parameters 
n=100000;                       %Number of simulation
sub_size = 125;
dt=gdouble(1/250);                       %Time step (1year = 250days)

%% Define storages
SimulPrices=repmat(CurrentPrice,n,1);
CorrWiener = gzeros(T-1,2,sub_size,'double');

%% Generating the paths of stock prices by Geometric Brownian Motion
UpperTriangle=gdouble(chol(Corr));    %UpperTriangle Matrix by Cholesky decomposition.  I don't have access to double precision Chol in Jacket so use CPU chol and convert.
Volr = repmat(Vol,[T-1,1,sub_size]);
Divr = repmat(Div,[T-1,1,sub_size]);

for ii = 1:sub_size:n
for i=1:sub_size
     CorrWiener(:,:,i)=grandn(T-1,2)*UpperTriangle;
end

%% do simulation
sim = cumprod(exp((r-Divr-Volr.^2./2).*dt+Volr.*sqrt(dt).*CorrWiener));
%get just the final prices
SimulPrices(ii:ii+sub_size-1,:) = SimulPrices(ii:ii+sub_size-1,:).*reshape(sim(end,:,:),2,sub_size)';
end
%% Plot the distribution of final prices
% Comment this section out if doing timings
%subplot(1,2,1);hist(SimulPrices(:,1),100);
%subplot(1,2,2);hist(SimulPrices(:,2),100);
gsync;