
close all
clear all

%
% Q1(a)
%

disp(' '); disp('Q1(a)'); disp(' ')

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

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

%
% Q1(b)
%

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

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(c) -- extra piece doing PCA factorisation
%

%disp(' '); disp('Q1(c)'); disp(' ')
%
%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(a)
%

disp(' '); disp('Q2(a)'); disp(' ')

M = 1000;  % number of independent Monte Carlo averages
N = 1000;  % number of samples in each Monte Carlo average

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

W = sqrt(T)*randn(N,M);
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);

%
% Q2(b)
%

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

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','NorthWest')

%
% Q3(a)
%

disp(' '); disp('Q3(a)'); disp(' ')

N = 10^6;

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

W  = sqrt(T)*randn(N,1);
fp = exp(-r*T)*max(0,S0*exp((r-0.5*sigma^2)*T+sigma*W)-K);
fm = exp(-r*T)*max(0,S0*exp((r-0.5*sigma^2)*T-sigma*W)-K);

fbar = sum(fp+fm)/(2*N);
sig2 = sum(fp.^2+fm.^2)/(2*N) - fbar^2;
cov  = sum((fp-fbar).*(fm-fbar))/N;
corr = cov/sig2;

variance_reduction_factor = 0.5*(1+corr);

fprintf('correlation:       %f\n',corr)
fprintf('variance reduction: %f\n',variance_reduction_factor)


%
% Q3(b)
%

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

N = 10^6;

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

W = sqrt(T)*randn(N,1);
S = S0*exp((r-0.5*sigma^2)*T+sigma*W);
f = exp(-r*T)*max(0,S-K);
g = exp(-r*T)*S;

fbar  = sum(f)/N;
gbar  = sum(g)/N;
sig2f = sum(f.^2)/N - fbar^2;
sig2g = sum(g.^2)/N - gbar^2;
cov   = sum(f.*g)/N - fbar*gbar;
corr  = cov / sqrt(sig2f*sig2g);

variance_reduction_factor = 1-corr^2;

fprintf('correlation:        %f\n',corr)
fprintf('variance reduction: %f\n',variance_reduction_factor)

%
% Q4(a)
%

disp(' '); disp('Q4(a)'); disp(' ')

N = 10^6;

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

W = sqrt(T)*randn(N,1);
S = S0*exp((r-0.5*sigma^2)*T+sigma*W);
f = exp(-r*T)*0.5*(1+sign(K-S));

fbar = sum(f)/N;
sig2 = sum(f.^2)/N - fbar^2;

fexact = exp(-r*T)*normcdf((log(K)-(log(S0)+(r-0.5*sigma^2)*T))/(sigma*sqrt(T)));

N_required = 9*sig2 / (0.1*fbar)^2;

fprintf('numerical value: %f  +/- %f\n',fbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',fexact)
fprintf('number of paths required: %g\n',N_required)


%
% Q4(b)
%

disp(' '); disp('Q4(b)'); disp(' ')

N = 10^6;

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

W     = sqrt(T)*randn(N,1);
mu1   = (r-0.5*sigma^2)*T;
mu2   = log(K) - log(S0);
shift = mu2 - mu1;
R     = exp(-(0.5*shift+sigma*W).*shift/(sigma^2*T));
S     = K*exp(sigma*W);
f     = exp(-r*T)*0.5*(1+sign(K-S)) .* R;

fbar = sum(f)/N;
sig2 = sum(f.^2)/N - fbar^2;

fexact = digital_put(r,sigma,T,S0,K,'value');

N_required = 9*sig2 / (0.1*fbar)^2;

fprintf('numerical value: %f  +/- %f\n',fbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',fexact)
fprintf('number of paths required: %g\n',N_required)

%
% Q5(a)
%

disp(' '); disp('Q5(a)')

N = 10^6;

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

W = sqrt(T)*randn(N,1);
S = S0*exp((r-0.5*sigma^2)*T+sigma*W);
f = exp(-r*T)*max(0,S-K);

%
% delta
%

g = f .* W/(S0*sigma*T);

gbar = sum(g)/N;
sig2 = sum(g.^2)/N - gbar^2;

gexact = european_call(r,sigma,T,S0,K,'delta');

disp(' '); disp('delta'); disp(' ')
fprintf('numerical value: %f  +/- %f\n',gbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',gexact)

%
% vega
%

g = f .* ( (W.^2-T)/(sigma*T) - W);

gbar = sum(g)/N;
sig2 = sum(g.^2)/N - gbar^2;

gexact = european_call(r,sigma,T,S0,K,'vega');

disp(' '); disp('vega'); disp(' ')
fprintf('numerical value: %f  +/- %f\n',gbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',gexact)

%
% Q5(b)
%

disp(' '); disp('Q5(b)')

fp = exp(-r*T)*0.5*(1+sign(S-K));

%
% delta
%

g = fp .* S/S0;

gbar = sum(g)/N;
sig2 = sum(g.^2)/N - gbar^2;

gexact = european_call(r,sigma,T,S0,K,'delta');

disp(' '); disp('delta'); disp(' ')
fprintf('numerical value: %f  +/- %f\n',gbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',gexact)

%
% vega
%

g = fp .* S.*(W-sigma*T);

gbar = sum(g)/N;
sig2 = sum(g.^2)/N - gbar^2;

gexact = european_call(r,sigma,T,S0,K,'vega');

disp(' '); disp('vega'); disp(' ')
fprintf('numerical value: %f  +/- %f\n',gbar,3*sqrt(sig2/N))
fprintf('analytic value:  %f\n',gexact)

