%
% Practical 3 -- Q3
%
% Test weak convergence of Euler method
% for European call options
%
% 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   = [80 90 100 110]';

Ve  = european_call(r,sig,T,S0,K,'value');

%
% Monte Carlo simulation comparing to exact solution
%

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

for p = 1:5
  N = 2^p;
  h = T/N;

  sum1 = 0;
  sum2 = 0;

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

    S = S0*ones(1,m2);

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

    P = exp(-r*T)*max(ones(size(K))*S-K*ones(size(S)),0);

    sum1 = sum1 + sum(P,2);
    sum2 = sum2 + sum(P.^2,2);
  end

  V  = sum1/M;
  sd = sqrt((sum2/M - (sum1/M).^2)/(M-1));

  hh(p)     = h;
  err1(:,p) = V-Ve;
  err2(:,p) = 3*sd;
end

figure;
loglog(hh,abs(err1(1,:)),'b-*',hh,err2(1,:),'r-*', ...
       hh,abs(err1(2,:)),'b-o',hh,err2(2,:),'r-o', ...
       hh,abs(err1(3,:)),'b-x',hh,err2(3,:),'r-x', ...
       hh,abs(err1(4,:)),'b-s',hh,err2(4,:),'r-s')
title('Weak convergence -- comparison to exact solution')
xlabel('h'); ylabel('Error');
legend(' K=80 error',' MC error', ...
       ' K=90 error',' MC error', ...
       ' K=100 error',' MC error', ...
       ' K=110 error',' MC error', 'Location','EastOutside')
axis([.02 0.5 0.005 2])

%
% Monte Carlo simulation comparing to 2h simulation
%

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

for p = 1:7
  N = 2^p;
  h = T/N;

  sum1 = 0;
  sum2 = 0;

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

    S  = S0*ones(1,m2);
    S2 = S0*ones(1,m2);

    for n = 1:N/2
      dW1 = sqrt(h)*randn(1,m2);
      S  = S.*(1+r*h+sig*dW1);    
      dW2 = sqrt(h)*randn(1,m2);
      S  = S.*(1+r*h+sig*dW2);

      S2 = S2.*(1+r*2*h+sig*(dW1+dW2));
    end

    P  = exp(-r*T)*max(ones(size(K))*S-K*ones(size(S)),0);
    P2 = exp(-r*T)*max(ones(size(K))*S2-K*ones(size(S)),0);

    sum1 = sum1 + sum(P-P2,2);
    sum2 = sum2 + sum((P-P2).^2,2);
  end

  V  = sum1/M;
  sd = sqrt((sum2/M - (sum1/M).^2)/(M-1));

  hh(p)     = h;
  err1(:,p) = V;
  err2(:,p) = 3*sd;
end

figure;
loglog(hh,abs(err1(1,:)),'b-*',hh,err2(1,:),'r-*', ...
       hh,abs(err1(2,:)),'b-o',hh,err2(2,:),'r-o', ...
       hh,abs(err1(3,:)),'b-x',hh,err2(3,:),'r-x', ...
       hh,abs(err1(4,:)),'b-s',hh,err2(4,:),'r-s')
title('Weak convergence -- difference from 2h approximation')
xlabel('h'); ylabel('Error');
legend(' K=80 error',' MC error', ...
       ' K=90 error',' MC error', ...
       ' K=100 error',' MC error', ...
       ' K=110 error',' MC error', 'Location','EastOutside')
axis([.005 0.5 0.001 0.5])
