
close all
clear all

%
% Q1(a)
%

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


%
% Q1(b)
%

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

%
% Q2(a)
%

disp(' '); disp('Q2(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;

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


%
% Q2(b)
%

disp(' '); disp('Q2(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 = 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;

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

%
% Q3
%

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

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

for option = 1:2
  if option==1
    gexact = european_call(r,sigma,T,S0,K,'delta');
    bumps  = 10.^linspace(-1,1,20);
  else
    gexact = european_call(r,sigma,T,S0,K,'vega');
    bumps  = 10.^linspace(-3,-1,20);
  end
  
for pass = 1:2
  vals  = [];
  vars  = [];

  for bump = bumps
    if (pass==2)
      rng('default');          % reset random number generator
    end

    W = sqrt(T)*randn(N,1);
    if (option==1)
      S = (S0-bump)*exp((r-0.5*sigma^2)*T+sigma*W);
    else
      S = S0*exp((r-0.5*(sigma-bump)^2)*T+(sigma-bump)*W);
    end
    fm = exp(-r*T)*max(0,S-K);

    if (pass==2)
      rng('default');          % reset random number generator
    end

    W = sqrt(T)*randn(N,1);
    if (option==1)
      S = (S0+bump)*exp((r-0.5*sigma^2)*T+sigma*W);
    else
      S = S0*exp((r-0.5*(sigma+bump)^2)*T+(sigma+bump)*W);
    end
    fp = exp(-r*T)*max(0,S-K);

    Del = (fp-fm)/(2*bump);

    val =  sum(Del)/N;
    var = (sum(Del.^2)/N - val^2)/(N-1);

    vals  = [vals  val];
    vars  = [vars  var];
  end

  figure(2*(option-1)+pass);
  pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);

  subplot(2,1,1)
  plot(bumps,vals-gexact,'-x');
  xlabel('bump');
  if option==1
    ylabel('Delta error')
  else
    ylabel('Vega error')
  end
  if (pass==1)
    title('different random numbers')
  else
    title('same random numbers')
  end

  subplot(2,1,2)
  plot(bumps,sqrt(vars),'-x');
  xlabel('bump'); ylabel('s.d.')

  if (pass==1)
    print('-deps2c','lec7_1.eps')
  else
    print('-deps2c','lec7_2.eps')
  end
end

end

%
% Q3 -- extra (LRM)
%

disp(' '); disp('Q3 -- extra (LRM)')

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(' ')
disp(sprintf('numerical value: %f  +/- %f',gbar,3*sqrt(sig2/N)))
disp(sprintf('analytic value:  %f',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(' ')
disp(sprintf('numerical value: %f  +/- %f',gbar,3*sqrt(sig2/N)))
disp(sprintf('analytic value:  %f',gexact))

%
% Q3 -- extra (pathwise sensitivity)
%

disp(' '); disp('Q3 -- extra (pathwise)')

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(' ')
disp(sprintf('numerical value: %f  +/- %f',gbar,3*sqrt(sig2/N)))
disp(sprintf('analytic value:  %f',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(' ')
disp(sprintf('numerical value: %f  +/- %f',gbar,3*sqrt(sig2/N)))
disp(sprintf('analytic value:  %f',gexact))

