close all; clear all;
rng('default')  % initialises random number generator

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

N = 1:1000000;
U = rand(1,max(N));    % uniform random variable
Y = ncfinv(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);

sum1 = cumsum(F);      % cumulative summation of
sum2 = cumsum(F.^2);   % payoff and its square
val  = sum1./N;
rms  = sqrt(sum2./N - val.^2);

err  = european_call(r,sig,T,S0,K,'value') - val;

plot(N,err,                ...
     N,err-3*rms./sqrt(N), ...
     N,err+3*rms./sqrt(N))
axis([0 length(N) -1 1])
xlabel('N'); ylabel('Error')
legend('MC error','lower bound','upper bound')

print('-depsc2','call_convergence.eps')

