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

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

%
% first figure
%

U1 = rand(1,50);        % uniform random variable
Y1 = norminv(U1);       % inverts Normal cum. fn.
S1 = S0*exp((r-sig^2/2)*T + sig*sqrt(T)*Y1);
F1 = exp(-r*T)*max(0,S1-K);

N  = 10000;
U2 = ((1:N)-0.5)/N;     % uniformly spaced values
Y2 = norminv(U2);       % inverts Normal cum. fn.
S2 = S0*exp((r-sig^2/2)*T + sig*sqrt(T)*Y2);
F2 = exp(-r*T)*max(0,S2-K);

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
plot(U1,F1,'rx',U2,F2,'b-')
title('Discounted payoff')
xlabel('U')

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

%
% second figure
%

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

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;

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*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','lec2b.eps')
