%
% Practical 4 - Question 3
% Use Euler-Maruyama method and Pathwise Sensitivity Method to
% estimate Vega for Down & Out call option based on usual GBM
% Use 'Alternative 1' to retain O(h) weak convergence
%
% Test problem:   dS   = r*S dt + sig*S dW
%

clear all; close all
randn('state',0)

%
% problem parameters and exact solution
%

r   = 0.05; 
sig = 0.5; 
T   = 1; 
S0  = 100; 
K   = 110;
B   = 90;

Vb_exact  = down_and_out_call(r,sig,T,S0,K,B,'value');
vegab_exact  = down_and_out_call(r,sig,T,S0,K,B,'vega');

%
% Monte Carlo simulation comparing to exact solution
%

M  = 1e+7;   % total number of Monte Carlo paths
M2 = 1e+4;   % number of paths at a time 
N = 64;

h = T/N;

sum1 = 0;
sum2 = 0;
pssum1 = 0;
pssum2 = 0;

for m = 1:M2:M
    m2 = min(M2,M-m+1); 

    S = S0*ones(1,m2);
    prob = ones(1,m2);
    s = zeros(1,m2);
    dPnsum = zeros(1,m2);

    for n = 1:N
      dW = sqrt(h)*randn(1,m2);
      sn = s;
      s  = s.*(1+r*h+sig*dW)+S.*dW;
      Sn = S;
      S  = S.*(1+r*h+sig*dW);
      SnmB = Sn-B; SmB = S-B;
      Pn = exp(-max(0,2*(S-B).*(Sn-B)./(sig^2*h*Sn.^2)));
      dPndsig = -Pn.*(2*sig*(s.*SnmB+SmB.*sn).*Sn - ...
                      4*sig*sn.*SmB.*SnmB - ...
                      4*SmB.*SnmB.*Sn) ./ (sig^3*h*Sn.^3);
      prob = prob.*(1-Pn);
      dPnsum = dPnsum + dPndsig./(1-Pn);
    end

    P    = exp(-r*T)*max(S-K,0).*prob;
    sum1 = sum1 + sum(P);
    sum2 = sum2 + sum(P.^2);

    dPnsum(prob==0)=0;
    psP    = exp(-r*T)*s.*(S>=K).*prob - P.*dPnsum;
    pssum1 = pssum1 + sum(psP);
    pssum2 = pssum2 + sum(psP.^2);
end

Vb_approx  = sum1/M;
sig2 = (sum2/M - (sum1/M)^2);

vegab_approx  = pssum1/M;
pssig2 = (pssum2/M - (pssum1/M)^2);

disp(sprintf('numerical value: %f  +/- %f',Vb_approx,3*sqrt(sig2/M)))
disp(sprintf('analytic value:  %f',Vb_exact))
disp(sprintf('numerical vega: %f  +/- %f',vegab_approx,3*sqrt(pssig2/M)))
disp(sprintf('analytic vega:  %f',vegab_exact))

