%
% 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;


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

% two passes for Cholesky and PCA factorisations of correlation matrix

  for opass = 1:2  

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

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

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

    for ipass = 1:3
      N  = 64;        % number of timesteps
      M2 = 32;        % number of randomisations
      M  = 2^15 / M2; % number of paths in each "family"

      h  = T/N;

      sum1 = 0;
      sum2 = 0;

      if ipass==1
        disp('Sobol with BB')
      elseif ipass==2
        disp('Sobol without BB')
      else
        disp('plain MC')
      end

      for m = 1:M2
        if ipass==1
%
%         Sobol points with Brownian Bridge construction of Brownian increments
%
          Ps = sobolset(5*N);    % dimension 5
          Ps = scramble(Ps,'MatousekAffineOwen');
          U  = net(Ps,M)';     % 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*N);    % dimension 5
          Ps = scramble(Ps,'MatousekAffineOwen');
          U  = net(Ps,M)';     % generate set of M Sobol points
          Z  = norminv(U);         % inverts Normal cum. fn.
          Z  = reshape(Z,5,N*M);
          dW = sqrt(h)*L*Z;
          dW = reshape(dW,5*N,M);

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

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

        for n = 1:N
          S  = S.*(1+r*h+sigma*dW(5*n-4:5*n,:));
          Save = Save + 0.2*sum(S,1)/N;
          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)/M;

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

      V  = sum1/M2;
      sd = sqrt((sum2/M2 - V^2)/(M2-1));
      if ipass==3
        fprintf(' MC_val      = %f \n',V)
        fprintf(' MC_std_dev  = %f \n\n',sd)
      else
        fprintf(' QMC_val      = %f \n',V)
        fprintf(' QMC_std_dev  = %f \n\n',sd)
      end
    end
  end
end
