S0=[50;55;35;45;5350];
X0=[1;1;1.49;1.28;1/101];
rUSD=0.004;rUK=0.004;rGermany=0.001;rJPY=0.001; % interest rates
y=0.02*ones(5,1); % dividend yields
T=2;
steps=256;

dt=T/steps;
sims=10000;
subSize = 250;
r=zeros(5,1);
r(1)=rUSD-y(1)-v(1)^2/2;
r(2)=rUSD-y(2)-v(2)^2/2;
r(3)=rUK-y(3)-R(3,6)*v(6)*v(3)-v(3)^2/2;
r(4)=rGermany-y(4)-R(4,7)*v(7)*v(4)-v(4)^2/2;
r(5)=rJPY-y(5)-R(5,8)*v(8)*v(5)-v(5)^2/2;

time_indices = [16 32 64 128 256];

K0 = [50;55;35;45;5350];
cap = 30;
i = 1;

corrmat = R;
upperTriangle=chol(corrmat(1:length(S0),1:length(S0)));

numBatches = 10;
latinHypercubePrices = zeros(numBatches,1);
for numBatch = 1:numBatches
    Spath=LatinHypercube(upperTriangle,r,S0,v(1:length(S0)),sims,steps,T); % latin
    [~, latinHypercubePrices(numBatch,1), standard_error] = Payoff(Spath, K0, time_indices, cap, X0, rUSD, T);
end
disp('Latin Hypercube');
disp(mean(latinHypercubePrices));
disp(std(latinHypercubePrices)/sqrt(numBatches));

sobolPrices = zeros(numBatches, 1);
sobolsequence=sobolset(size(R,2)-2, 'Skip', 1e3);
for numBatch = 1:numBatches
    SobolPoints=sobolsequence((numBatch-1)*sims+1:numBatch*sims,1:size(R,2)-2);
    Spath = Sobol(upperTriangle,r,S0,v(1:length(S0)),sims,steps,T, SobolPoints);
    [~, sobolPrices(numBatch,1), standard_error] = Payoff(Spath, K0, time_indices, cap, X0, rUSD, T);
end
disp('Sobol');
disp(mean(sobolPrices));
disp(std(sobolPrices)/sqrt(numBatches));

momentMatchingPrices = zeros(numBatches, 1);
for numBatch = 1:numBatches
    Spath=MomentMatching(upperTriangle,r,S0,v(1:length(S0)),sims,steps,T, subSize); % Moment matching
    [~, momentMatchingPrices(numBatch, 1), standard_error] = Payoff(Spath, K0, time_indices, cap, X0, rUSD, T);
end
disp('MomentMatching')
disp(mean(momentMatchingPrices));
disp(std(momentMatchingPrices)/sqrt(numBatches));

Spath=SamplePath2(upperTriangle,r,S0,v(1:length(S0)),sims,steps,T, subSize); % naive MC
[payoffs_Naive, naivePrice, standard_error] = Payoff(Spath, K0, time_indices, cap, X0, rUSD, T);
disp('Naive');
disp(naivePrice);
disp(standard_error);


[Spath1,Spath2]=AntitheticVar(upperTriangle,r,S0,div,v(1:length(S0)),sims,steps,T,subSize);
[antitheticPrice, standard_error] = AntitheticPayoff(Spath1,Spath2, K0, time_indices, cap, X0, rUSD, T);
disp('Antithetic');
disp(antitheticPrice);
disp(standard_error);