function [SimulPrices] = GPU_PCT_corr3( n,sub_size)
%GPU_PCT_corr3 - A third attempt at a GPU version of a simulation of two correlated assets
%n - number of simulations
%sub_size - number of simulations per chunk
%% Correlated asset information
CurrentPrice = gpuArray([78 102]);       %Initial Prices of the two stocks
Corr = gpuArray([1 0.4; 0.4 1]);         %Correlation Matrix
T = gpuArray(500);                       %Number of days to simulate = 2years = 500days 
Div=gpuArray([0.01 0.01]);               %Dividend
Vol=gpuArray([0.2 0.3]);                 %Volatility

%% Market Information
r = 0.03;                      %Risk-free rate

%% Simulation parameters 
dt=1/250;                       %Time step (1year = 250days)

%% Define storages
SimulPrices=repmat(CurrentPrice,n,1);
CorrWiener = parallel.gpu.GPUArray.zeros([T-1,2,sub_size]);

%% Generating the paths of stock prices by Geometric Brownian Motion
UpperTriangle=chol(Corr);    %UpperTriangle Matrix by Cholesky decomposition
Volr = repmat(Vol,[T-1,1,sub_size]);
Divr = repmat(Div,[T-1,1,sub_size]);

CorrWiener_final = parallel.gpu.GPUArray.zeros(T-1,2,sub_size);
for ii = 1:sub_size:n
%Generate correlated random numbers
randoms = parallel.gpu.GPUArray.randn(sub_size*(T-1),2);
CorrWiener = randoms*UpperTriangle;

%Put numbers into n pages of T-1*stocks
%Each page represents a separate simulation
CorrWiener = reshape(CorrWiener,(T-1),sub_size,2);
%poor man's permute since GPU permute if not available in 2011b
for s = 1:2
    CorrWiener_final(:, s, :) = CorrWiener(:, :, s);
end

%% do simulation
sim = cumprod(exp((r-Divr-Volr.^2./2).*dt+Volr.*sqrt(dt).*CorrWiener_final));
%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);
end