%
% Practical 4 - Question 2
% Use Euler-Maruyama method and Pathwise Sensitivity Method to
% estimate Vega for European call option based on usual GBM
%
% Test problem:   dS   = r*S dt + sig*S dW
%

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

%
% problem parameters and exact solution
%

eps = 1e-20;
r   = 0.05; 
sig = 0.5 + 1i*eps; 
T   = 1; 
S0  = 100; 
K   = 110;

V_exact  = european_call(r,sig,T,S0,K,'value');
vega_exact  = european_call(r,sig,T,S0,K,'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 = 256;

h = T/N;

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

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

    S = S0*ones(1,m2);
    s = zeros(1,m2);
    
    for n = 1:N
      dW = sqrt(h)*randn(1,m2);
      s  = s.*(1+r*h+sig*dW)+S.*dW;    
      S  = S.*(1+r*h+sig*dW);    
    end

    % Check our error in dS_N\dsig
    %
    % s will have gained an imag part because of our CVT - we only want
    % real part to give us dS_N\dsig. 
    %
    % imag(S)/eps will give us CVT value for dS_N\dsig
    s = real(s); 
    err = max(err,max(abs(s - imag(S)/eps)));
    
    % remove imaginary part of S - no longer needed
    S = real(S);
    P = exp(-r*T)*max(S-K,0);
    psP = exp(-r*T)*(S>=K).*s;

    sum1 = sum1 + sum(P);
    sum2 = sum2 + sum(P.^2);
    pssum1 = pssum1 + sum(psP);
    pssum2 = pssum2 + sum(psP.^2);
end

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

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

disp(sprintf('numerical value: %f  +/- %f',V_approx,3*sqrt(sig2/M)))
disp(sprintf('analytic value:  %f',V_exact))
disp(sprintf('numerical vega: %f  +/- %f',vega_approx,3*sqrt(pssig2/M)))
disp(sprintf('analytic vega:  %f',vega_exact))
disp(['error in dS_N\dsig vs CVT result: ' num2str(err)]);

