%
% Evaluate error for European call option using
% Monte Carlo, Latin Hypercube and Sobol QMC sampling
% with either Cholesky or PCA factorisation
%

clear all; close all

rng('default')

for pass = 1:2  % two passes for Cholesky and PCA factorisations

  N = 2^20;

  r     = 0.05;
  sigma = 0.2;
  T     = 1;
  K     = 100;
  S0    = 100;

  Sigma = sigma^2*T*( eye(5) + 0.1*(ones(5)-eye(5)));

  if pass==1
    disp('Cholesky factorisation:'); disp(' ');
    L = chol(Sigma,'lower');
  else
    disp('PCA factorisation:'); disp(' ');
    [V,D] = eig(Sigma);
    L = V*diag(sqrt(diag(D)));
  end

% first do standard MC

  U = rand(5,N);         % uniform random variable
  Y = norminv(U);        % inverts Normal cum. fn.
  S = S0*exp((r-sigma^2/2)*T + L*Y);
  F = exp(-r*T)*max(0,sum(S,1)/5-K);

  sum1 = sum(F);
  sum2 = sum(F.^2);
  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/N);
  fprintf(' MC_val        = %f \n',val)
  fprintf(' MC_err_bound  = %f \n\n',err_bound)

% now do Latin Hypercube

  M = 1024;  % number of strata
  N = 1024;  % number of Latin Hypercube samples

  sum1 = 0;
  sum2 = 0;
  m = zeros(5,M);

  for n = 1:N
    for d = 1:5
      m(d,:) = randperm(M); % Latin Hypercube locations
    end

    U = (m-1+rand(5,M))/M;  % random location within each cube 
    Y = norminv(U);         % inverts Normal cum. fn.
    S = S0*exp((r-sigma^2/2)*T + L*Y);
    F = exp(-r*T)*max(0,sum(S,1)/5-K);

    sum1 = sum1 +  sum(F)/M;
    sum2 = sum2 + (sum(F)/M)^2;
  end

  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/N);
  fprintf(' LH_val        = %f \n',val)
  fprintf(' LH_err_bound  = %f \n\n',err_bound)

% now do QMC

  Ntot = 2^20;      % total number of samples to take
  N    = 2^4;       % number of sets
  M    = Ntot / N;  % number of QMC points
  
  sum1 = 0;
  sum2 = 0;
  
  for n = 1:N            % simulate for each set of Sobol points
    Ps = sobolset(5);    % dimension 5
    Ps = scramble(Ps,'MatousekAffineOwen');
    U  = net(Ps,M)';     % generate set of M Sobol points

    Y = norminv(U);      % inverts Normal cum. fn.
    S = S0*exp((r-sigma^2/2)*T + L*Y);
    F = exp(-r*T)*max(0,sum(S,1)/5-K);

    sum1 = sum1 +  sum(F)/M;
    sum2 = sum2 + (sum(F)/M)^2;
  end

  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/(N-1));
  fprintf(' QMC_val       = %f \n',val)
  fprintf(' QMC_err_bound = %f \n\n\n',err_bound)
end