%
% Evaluate error for European call option using Sobol points
%

clear all; close all

rng('default')

%
% 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:2

  for p = 2:12
    M = 2^p;  % number of points in each "family"
    N = 64;   % number of "families"

    for n = 1:N

% standard Monte Carlo
      if pass==1
        U  = rand(M,1);

% Sobol QMC
      else
        Ps = sobolset(1);     % dimension 1
        Ps = scramble(Ps,'MatousekAffineOwen');
        U  = net(Ps,M)';
      end

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

      Fave(n) = sum(F)/M;
    end
 
    V  = sum(Fave)/N;
    sd = sqrt((sum(Fave.^2)/N - (sum(Fave)/N)^2)/(N-1));

    Ms(p)   = M*N;
    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-*')
  if pass==1
    xlabel('N'); ylabel('Error');
    legend(' Error',' MC error bound','Location','NorthEast')
    title('MC convergence')
    print('-deps2c','lec5c.eps')
  else
    xlabel('M*N'); ylabel('Error');
    legend(' Error',' QMC error bound','Location','NorthEast')
    title('QMC convergence')
    print('-deps2c','lec5d.eps')
  end

end
