
close all
clear all

%
% Q1(b)
%

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

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

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

%
% Q1(c)
%

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

Sig = [ 4 1 1; 1 4 1; 1 1 4];

L = chol(Sig,'lower');
X = randn(3,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; 1 4 1; 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)*ncfinv(C),C)
xlabel('Y'); ylabel('C')
legend('Monte Carlo results','Central Limit Theorem',4)


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

%
% 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)
xlabel('N'); ylabel('value')
legend('mean','lower','upper','exact')

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

%
% 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);

disp(sprintf('correlation:       %f',corr))
disp(sprintf('variance reduction: %f',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;

disp(sprintf('correlation:        %f',corr))
disp(sprintf('variance reduction: %f',variance_reduction_factor))

%
% Q4(a)
%

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

for pass = 1:2  % two passes for Cholesky and PCA factorisations

  N = 10^6;

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

  Sigma = sigma^2*T*( eye(5) + 0.1*(ones(5)-eye(5)));

  if pass==1
    disp('Cholesky factorisation:'); disp(' ');
    L = chol(Sigma,'lower');
  else
    disp('PCA factorisation:'); disp(' ');
    [V,D] = eig(Sigma);
    L = V*diag(sqrt(diag(D)));
  end

% first do standard MC

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

  sum1 = sum(F);
  sum2 = sum(F.^2);
  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/N);
  fprintf(' MC_val        = %f \n',val)
  fprintf(' MC_err_bound  = %f \n\n',err_bound)

% now do Latin Hypercube

  M = 1000;  % number of strata
  N = 1000;  % number of Latin Hypercube samples

  sum1 = 0;
  sum2 = 0;
  m = zeros(5,M);

  for n = 1:N
    for d = 1:5
      m(d,:) = randperm(M); % Latin Hypercube locations
    end

    U = (m-1+rand(5,M))/M;  % random location within each cube 
    Y = ncfinv(U);          % inverts Normal cum. fn.
    S = S0*exp((r-sigma^2/2)*T + L*Y);
    F = exp(-r*T)*max(0,sum(S,1)/5-K);

    sum1 = sum1 +  sum(F)/M;
    sum2 = sum2 + (sum(F)/M)^2;
  end

  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/N);
  fprintf(' LH_val        = %f \n',val)
  fprintf(' LH_err_bound  = %f \n\n',err_bound)

% now do QMC -- code originally due to Dirk Nuyens

  Ntot = pow2(20);  % total number of samples to take
  N    = pow2(4);   % number of shifts
  M    = Ntot / N;  % number of QMC points to use
  
  shifts = rand(5, N);  % the 5-dimensional random shifts
  latticeseq_b2 init0;  % initialize QMC generator
  QMCsamples = latticeseq_b2(5, M);  % QMC samples

  sum1 = 0;
  sum2 = 0;
  
  for n = 1:N          % now simulate for each of the shifts
    U = mod1shift(QMCsamples, shifts(:,n)); % shift the QMC samples
    Y = ncfinv(U);         % inverts Normal cum. fn.
    S = S0*exp((r-sigma^2/2)*T + L*Y);
    F = exp(-r*T)*max(0,sum(S,1)/5-K);

    sum1 = sum1 +  sum(F)/M;
    sum2 = sum2 + (sum(F)/M)^2;
  end

  val  = sum1/N;
  var  = sum2/N - val^2;
  err_bound = 3*sqrt(var/N);
  fprintf(' QMC_val       = %f \n',val)
  fprintf(' QMC_err_bound = %f \n\n\n',err_bound)
end
