%
% Evaluate error for European call and barrier options
% using Monte Carlo and Sobol QMC sampling with either
% Cholesky or Brownian Bridge factorisation in time
%

clear all; close all

rng('default')

r      = 0.02;
sigma  = 0.2;
xi     = 0.2;
lambda = 2;
rho    = -0.2;
T      = 1;
K      = 100;
B      = 80;
S0     = 100;
v0     = 0.04;

for option = 1:2
  disp('-----------------')
  if option==1
    disp('Call option');
  elseif option==2
    disp('Barrier option');
  end
  disp('-----------------')

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

% Cholesky factorisation of correlation matrix
  L = chol(Sigma,'lower');
% alternative PCA factorisation of correlation matrix
%  [V,D] = eig(Sigma);
%  L = V*diag(sqrt(diag(D)));

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

  for pass = 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 pass==1
      disp('Sobol with BB')
    elseif pass==2
      disp('Sobol without BB')
    else
      disp('plain MC')
    end

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

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

      S    = repmat(S0,1,M);
      v    = repmat(v0,1,M);
      prob = repmat(1,1,M);

      for n = 1:N
	sig  = sqrt(abs(v));
        Sn   = S;
        S    = S.*(1+r*h+sig.*dW(2*n-1,:));
        v    = v + lambda*(sigma^2-v)*h + xi*sig.*dW(2*n,:);
        prob = prob.*(1 - exp(-max(0,2*(S-B).*(Sn-B)./(h*sig.^2.*Sn.^2))));
       end

      if option==1
        P = exp(-r*T)*max(S-K,0);
      elseif option==2
        P = exp(-r*T)*max(S-K,0).*prob;
      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 pass==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
