%
% Practical 4 - Question 1
% Use Euler-Maruyama method and Likelihood Ratio 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
%

r   = 0.05; 
sig = 0.5; 
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 = 64;

h = T/N;

sum1 = 0;
sum2 = 0;
lrsum1 = 0;
lrsum2 = 0;

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

    S = S0*ones(1,m2);
    z2sum = 0;
    
    for n = 1:N
      Z = randn(1,m2);
      S  = S.*(1+r*h+sig*sqrt(h)*Z);    
      z2sum = z2sum + Z.^2;
    end

    P = exp(-r*T)*max(S-K,0);
    lrP = P.*(z2sum-N)/sig;

    sum1 = sum1 + sum(P);
    sum2 = sum2 + sum(P.^2);
    lrsum1 = lrsum1 + sum(lrP);
    lrsum2 = lrsum2 + sum(lrP.^2);
end

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

vega_approx  = lrsum1/M;
lrsig2 = (lrsum2/M - (lrsum1/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(lrsig2/M)))
disp(sprintf('analytic vega:  %f',vega_exact))

