
close all; clear all;
rng('default')

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

for M = [1 10 100 1000]
  L = N/M;  ave=0; var=0;

  for m = 1:M
    U = (m-1+rand(1,L))/M;   % uniform random variable
    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);
    ave1 =  sum(F)/L;
    var1 = (sum(F.^2)/L - ave1^2)/(L-1);
    ave = ave + ave1/M;
    var = var + var1/M^2;
  end

  err = european_call(r,sig,T,S0,K,'value') - ave;
  disp(sprintf('\nM = %d, L = %d',M, L))
  disp(sprintf('error       = %f',err))
  disp(sprintf('error bound = %f',3*sqrt(var)))
end
