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

clear all; close all

rng('default')

r     = 0.02;
sigma = 0.2;
rho   = 0.2;
T     = 1;
K     = 100;
S0    = 100;

Sigma = eye(5) + rho*(ones(5)-eye(5));

[V,D] = eig(Sigma);
D = fliplr(flipud(D));    % flip so eigenvalues in descending order
V = fliplr(V);            % corresponding flip of eigenvectors
L = V*diag(sqrt(diag(D)));


for option = 1:3
  if option==1
    disp('Basket call option');
  elseif option==2
    disp('Basket Asian call option');
  elseif option==3
    disp('Basket lookback call option');
  end

% two passes for Cholesky and PCA factorisations of correlation matrix

  for M = 1:7
    N  = 2^(M+6);   % number of paths in each "family"
    M2 = 32;        % number of randomisations
    NT = 64;        % number of timesteps
    h  = T/NT;

% three inner passes for Cholesky and Brownian Bridge factorisations
% of covariance matrix in time, and use of plain MC

    for ipass = 1:3

      sum1 = 0;
      sum2 = 0;

      for m = 1:M2
        if ipass==1
%
%         Sobol points with Brownian Bridge construction of Brownian increments
%
          Ps = sobolset(5*NT);    % dimension 5
          Ps = scramble(Ps,'MatousekAffineOwen');
          U  = net(Ps,N)';     % generate set of M Sobol points
          Z  = norminv(U);         % inverts Normal cum. fn.
          dW = bb(Z,T,L);

        elseif ipass==2
%
%         Sobol points without Brownian Bridge construction
%
          Ps = sobolset(5*NT);    % dimension 5
          Ps = scramble(Ps,'MatousekAffineOwen');
          U  = net(Ps,N)';     % generate set of M Sobol points
          Z  = norminv(U);         % inverts Normal cum. fn.
          Z  = reshape(Z,5,NT*N);
          dW = sqrt(h)*L*Z;
          dW = reshape(dW,5*NT,N);

        else
%
%       standard random number generation
%
          dW = sqrt(h)*L*randn(5,NT*N);
          dW = reshape(dW,5*NT,N);
        end

        S    = repmat(S0,5,N);
        Save = zeros(1,N);
        Smax = repmat(S0,1,N);

        for n = 1:NT
          S  = S.*(1+r*h+sigma*dW(5*n-4:5*n,:));
          Save = Save + 0.2*sum(S,1)/NT;
          Smax = max(Smax,0.2*sum(S,1));
        end

        S = 0.2*sum(S,1); 

        if option==1
          P = exp(-r*T)*max(S-K,0);
        elseif option==2
          P = exp(-r*T)*max(Save-K,0);
        else
          P = exp(-r*T)*max(Smax-K,0);
        end

        P = sum(P)/N;

        sum1 = sum1 + sum(P);
        sum2 = sum2 + sum(P.^2);
      end

      V  = sum1/M2;
      sd = sqrt((sum2/M2 - V^2)/(M2-1));

      sd_array(M,ipass) = sd;
      NN_array(M,ipass) = N;
    end

  end

  figure;
  loglog(NN_array,sd_array);
  xlabel('N'); ylabel('std dev');
  legend('QMC BB','QMC no BB','MC');

  if option==1
    title('Basket call option'); 
  elseif option==2
    title('Basket Asian call option');
  elseif option==3
    title('Basket lookback call option');
  end

end
