
close all
clear all

%
% Q1(b)
%

disp(' ')
disp('Q1(b)')

N = 10^6;
U = rand(1,N);
X = norminv(U);

mean = sum(X)/N
var  = X*X'/N - mean^2

%
% Q1(c)
%

disp(' ')
disp('Q1(c)')

Sig = [ 4 1; 1 4 ];

L = chol(Sig,'lower');
X = randn(2,N);
Y = L*X;

mean = sum(Y,2)/N
var  = Y*Y'/N - mean*mean'

%
% Q1(d)
%

disp(' ')
disp('Q1(d)')

Sig = [ 4 1; 1 4 ];

[V,D] = eig(Sig);
L = V*diag(sqrt(diag(D)));
Y = L*X;

mean = sum(Y,2)/N
var  = Y*Y'/N - mean*mean'

%
% Q2(b)
%

disp(' ')
disp('Q2(b)')

M = 1000;
N = 1000;

fbar = - 2/pi^2;
sig2 = 1/(4*pi^2) + 1/6 - fbar^2

U = rand(N,1000);
f = U.*cos(pi*U);

%
% Q2(c)
%

disp(' ')
disp('Q2(c)')

Y = sort(sum(f,1)/N);
C = ((1:M) - 0.5)/M;

figure
plot(Y,C, fbar+sqrt(sig2/N)*norminv(C),C)
xlabel('Y'); ylabel('C')
legend('Monte Carlo results','Central Limit Theorem','Location','SouthEast')


figure
plot(C,Y, C,fbar+sqrt(sig2/N)*norminv(C))
ylabel('Y'); xlabel('C')
legend('Monte Carlo results','Central Limit Theorem','Location','SouthEast')

%
% Q2(d)
%

disp(' ')
disp('Q2(d)')

N = 10^6;
U = rand(1,N);
f = U.*cos(pi*U);

S1 = cumsum(f);
S2 = cumsum(f.^2);

NN   = 10^3:10^3:10^6;
mean = S1(NN)./NN;
sd   = sqrt( (S2(NN)./NN - mean.^2)./(NN-1) );

figure
plot(NN,mean, NN,mean-3*sd, NN,mean+3*sd, NN,fbar*ones(size(NN)))
xlabel('N'); ylabel('value')
legend('mean','lower','upper','exact')

figure
plot(NN,NN.*sd.^2, NN,sig2*ones(size(NN)))
xlabel('N'); ylabel('variance')
legend('estimated','exact')


%
% Q3(b)
%

disp(' ')
disp('Q3(b)')

M = 1000;
N = 1000;

r     = 0.05;
sigma = 0.2;
T     = 1;
K     = 100;
S0    = 100;

U = rand(N,1000);
W = sqrt(T)*norminv(U);
f = exp(-r*T)*max(0,S0*exp((r-0.5*sigma^2)*T+sigma*W)-K);

fbar = european_call(r,sigma,T,S0,K,'value');
sig2 = sum(sum((f-fbar).^2)) / (M*N);

%
% Q3(c)
%

disp(' ')
disp('Q3(c)')

Y = sort(sum(f,1)/N);
C = ((1:M) - 0.5)/M;

figure
plot(Y,C, fbar+sqrt(sig2/N)*norminv(C),C)
xlabel('Y'); ylabel('C')
legend('Monte Carlo results','Central Limit Theorem','Location','SouthEast')


figure
plot(C,Y, C,fbar+sqrt(sig2/N)*norminv(C))
ylabel('Y'); xlabel('C')
legend('Monte Carlo results','Central Limit Theorem','Location','SouthEast')

%
% Q3(d)
%

disp(' ')
disp('Q3(d)')

N = 10^6;
U = rand(1,N);
W = sqrt(T)*norminv(U);
f = exp(-r*T)*max(0,S0*exp((r-0.5*sigma^2)*T+sigma*W)-K);

S1 = cumsum(f);
S2 = cumsum(f.^2);

NN   = 10^3:10^3:10^6;
mean = S1(NN)./NN;
sd   = sqrt( (S2(NN)./NN - mean.^2)./(NN-1) );

figure
plot(NN,mean, NN,mean-3*sd, NN,mean+3*sd, NN,fbar*ones(size(NN)))
xlabel('N'); ylabel('value')
legend('mean','lower','upper','exact')

