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

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

Sigma = sig^2*T*( eye(5) + 0.1*(ones(5)-eye(5)));
L     = chol(Sigma,'lower');

N = 1:1000000;
U = rand(5,max(N));    % uniform random variable
Y = ncfinv(U);         % inverts Normal cum. fn.
S = S0*exp((r-sig^2/2)*T + L*Y);
F = exp(-r*T)*max(0,sum(S,1)/5-K);

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

plot(N,val,                ...
     N,val-3*rms./sqrt(N), ...
     N,val+3*rms./sqrt(N))
axis([0 length(N) 14.5 16])
xlabel('N'); ylabel('Value')
legend('MC value','lower bound','upper bound')

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

