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

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

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

%
% first figure
%

N = 1:10^6;
Y = randn(1,max(N));    % Normal random variable
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;
sd  = sqrt(sum2./N - val.^2);

err = val - exact;

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
plot(N,err, N,err-3*sd./sqrt(N), N,err+3*sd./sqrt(N))
axis([0 length(N) -1 1])
xlabel('N'); ylabel('Error')
legend('MC error','lower bound','upper bound')

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

%
% second set of CLT figures
%

for pass = 1:3

if pass==1
  M = 100;   % number of independent tests
  N = 100;   % number of samples within each test
elseif pass==2
  M = 10000;
  N = 100;
elseif pass==3
  M = 10000;
  N = 10000;
end

Y = randn(M,N);    % Normal random variable
S = S0*exp((r-sig^2/2)*T + sig*sqrt(T)*Y);
F = exp(-r*T)*max(0,S-K);

sig2 = sum(sum((F-exact).^2)) / (M*N);

F = sum(F,2) / N;
F = sort([F; F]);  % repeat each entry
C = floor(0.5*(1:2*M))' / M;

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
plot(F,C, exact+sqrt(sig2/N)*norminv(C),C)
xlabel('MC value'); ylabel('CDF')
legend('MC results','CLT','Location','SouthEast')

if pass<3
  axis([10 25 0 1]);
end

print('-depsc2',sprintf('clt_%1d.eps',pass))

end
