%
% Test weak convergence of Euler method for 
% European call option using Sobol points
%
% Test problem:   dS   = r*S dt + sig*S dW
%

clear all; close all

%
% problem parameters and exact solution
%

r   = 0.05; 
sig = 0.5; 
T   = 1; 
S0  = 100; 
K   = 110;

Ve  = european_call(r,sig,T,S0,K,'value');

for pass = 1:3

%
% Quasi-Monte Carlo simulation comparing to exact solution
%

  N  = 128;   % number of timesteps
  M2 = 64;    % number of scrambles

  h  = T/N;

  for p = 2:12
    p
    M = 2^p;  % number of paths in each "family"

    sum1 = 0;
    sum2 = 0;

    for m = 1:M2
%
% Sobol quasi-random number generation with scrambling,
% and Brownian Bridge construction of Brownian increments
%
      if pass==1
        Ps = sobolset(N);
        Ps = scramble(Ps,'MatousekAffineOwen');

        U  = net(Ps,M)';
        Z  = norminv(U);
        dW = bb(Z,T);

%
% Sobol quasi-random number generation with scrambling,
% but no Brownian Bridge construction of Brownian increments
%
      elseif pass==2
        Ps = sobolset(N);
        Ps = scramble(Ps,'MatousekAffineOwen');

        U  = net(Ps,M)';
        Z  = norminv(U);
        dW = sqrt(h)*Z;

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

      S = repmat(S0,1,M);

      for n = 1:N
        S  = S.*(1+r*h+sig*dW(n,:));    
      end

      P = exp(-r*T)*max(S-K,0);
      P = sum(P)/M;

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

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

    Ms(p)   = M;
    err1(p) = V-Ve;
    err2(p) = 3*sd;
  end

  figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
  loglog(Ms,abs(err1),'b-*',Ms,err2,'r-*')
  xlabel('N'); ylabel('Error');
  legend(' Error',' MC error bound','Location','northeast')
  axis([1 4096 0.005 10])
  if pass==1
    title('comparison to exact solution')
    print('-deps2c','lec13a.eps')
  elseif pass==2
    title('comparison to exact solution')
    print('-deps2c','lec13b.eps')
  else
    title('comparison to exact solution')
    print('-deps2c','lec13c.eps')
  end

end
