%
% Test weak convergence of Euler method
% for barrier and lookback options
%
% Test problem: dS = r*S dt + sig*S dW
%

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

disp(' ')
disp('WARNING: this might take 5-15 mins to run!')
disp(' ')

%
% problem parameters and exact solution
%

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

Vb  = down_and_out_call(r,sig,T,S0,K,B,'value');
Vl  = lookback_call(r,sig,T,S0,'value');

%
% Monte Carlo simulation comparing to exact solution
%

if 1

% pass=1 no correction, pass=2 with correction

for pass = 1:2

if pass==1
  M  = 1e+6;   % total number of Monte Carlo paths
else
  M  = 1e+8;   % total number of Monte Carlo paths
end
M2 = 1e+4;   % number of paths at a time 

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

  sum1 = 0;
  sum2 = 0;
  sum3 = 0;
  sum4 = 0;

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

    S = S0*ones(1,m2);
    prob = ones(1,m2);
    Mn = S;

    for n = 1:N
      dW = sqrt(h)*randn(1,m2);
      Sn = S;
      S  = S.*(1+r*h+sig*dW);
      if pass==1
	Mn = min(Mn,S);
      else
      prob = prob.*(1 - exp(-max(0,2*(S-B).*(Sn-B)./(sig^2*h*Sn.^2))));
      Mn = min(Mn, 0.5*(S+Sn-sqrt( (S-Sn).^2-2*sig^2*h*Sn.^2.*log(rand(1,m2)))) );
      end
    end

    if pass==1
      P  = exp(-r*T)*max(S-K,0)*0.5.*(sign(Mn-B)+1);
    else
      P  = exp(-r*T)*max(S-K,0).*prob;
    end
    sum1 = sum1 + sum(P);
    sum2 = sum2 + sum(P.^2);

    P    = exp(-r*T)*(S-Mn);
    sum3 = sum3 + sum(P);
    sum4 = sum4 + sum(P.^2);
  end

  hh(p) = h;

  V  = sum1/M;
  sd = sqrt((sum2/M - (sum1/M)^2)/(M-1));
  err1(p) = V-Vb;
  err2(p) = 3*sd;

  V  = sum3/M;
  sd = sqrt((sum4/M - (sum3/M)^2)/(M-1));
  err3(p) = V-Vl;
  err4(p) = 3*sd;
end

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
loglog(hh,abs(err1),'b-*',hh,err2,'r-*')
%title('Barrier -- comparison to exact solution')
xlabel('h'); ylabel('Error');
legend(' Weak error',' MC error',2)

if pass == 1
  axis([.01 0.5 0.01 50])
  print('-deps2c','barrier_1.eps')
else
  axis([.01 0.5 0.001 5])
  print('-deps2c','barrier_3.eps')
end

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
loglog(hh,abs(err3),'b-*',hh,err4,'r-*')
%title('Lookback -- comparison to exact solution')
xlabel('h'); ylabel('Error');
legend(' Weak error',' MC error',2)

if pass==1
  axis([.01 0.5 0.05 50])
  print('-deps2c','lookback_1.eps')
else
  axis([.01 0.5 0.005 20])
  print('-deps2c','lookback_3.eps')
end

end

end

%
% Monte Carlo simulation comparing to 2h simulation
%

if 1

for pass = 1:2

if pass==1
  M  = 1e+5;   % total number of Monte Carlo paths
else
  M  = 1e+7;   % total number of Monte Carlo paths
end
M2 = 1e+4;   % number of paths at a time 

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

  sum1 = 0;
  sum2 = 0;
  sum3 = 0;
  sum4 = 0;

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

    S  = S0*ones(1,m2);
    S2 = S0*ones(1,m2);
    prob  = ones(1,m2);
    prob2 = ones(1,m2);
    Mn  = S;
    Mn2 = S2;

    for n = 1:N/2
      dW1 = sqrt(h)*randn(1,m2);
      Sn = S;
      S  = S.*(1+r*h+sig*dW1);
      if pass==1
        Mn  = min(Mn,S);
      else
      prob = prob.*(1 - exp(-max(0,2*(S-B).*(Sn-B)./(sig^2*h*Sn.^2))));
      Mn = min(Mn, 0.5*(S+Sn-sqrt( (S-Sn).^2-2*sig^2*h*Sn.^2.*log(rand(1,m2)))) );
      end

      dW2 = sqrt(h)*randn(1,m2);
      Sn = S;
      S  = S.*(1+r*h+sig*dW2);
      if pass==1
        Mn  = min(Mn,S);
      else
      prob = prob.*(1 - exp(-max(0,2*(S-B).*(Sn-B)./(sig^2*h*Sn.^2))));
      Mn = min(Mn, 0.5*(S+Sn-sqrt( (S-Sn).^2-2*sig^2*h*Sn.^2.*log(rand(1,m2)))) );
      end

      Sn2 = S2;
      S2 = S2.*(1+r*2*h+sig*(dW1+dW2));
      if pass==1
        Mn2  = min(Mn2,S2);
      else
      prob2 = prob2.*(1 - exp(-max(0,2*(S2-B).*(Sn2-B)./(sig^2*2*h*Sn2.^2))));
      Mn2 = min(Mn2, 0.5*(S2+Sn2-sqrt( (S2-Sn2).^2-2*sig^2*2*h*Sn2.^2.*log(rand(1,m2)))) );
      end
    end

    if pass==1
      P    = exp(-r*T)*max(S-K,0)*0.5.*(sign(Mn-B)+1);
      P2   = exp(-r*T)*max(S2-K,0)*0.5.*(sign(Mn2-B)+1);
    else
      P    = exp(-r*T)*max(S-K,0).*prob;
      P2   = exp(-r*T)*max(S2-K,0).*prob2;
    end
    sum1 = sum1 + sum(P-P2);
    sum2 = sum2 + sum((P-P2).^2);

    P    = exp(-r*T)*(S-Mn);
    P2   = exp(-r*T)*(S2-Mn2);
    sum3 = sum3 + sum(P-P2);
    sum4 = sum4 + sum((P-P2).^2);
  end

  hh(p) = h;

  Vd = sum1/M;
  sd = sqrt((sum2/M - (sum1/M)^2)/(M-1));
  err1(p) = Vd;
  err2(p) = 3*sd;

  Vd = sum3/M;
  sd = sqrt((sum4/M - (sum3/M)^2)/(M-1));
  err3(p) = Vd;
  err4(p) = 3*sd;
end

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
loglog(hh,abs(err1),'b-*',hh,err2,'r-*')
%title('Barrier -- comparison to 2h approximation')
xlabel('h'); ylabel('Error');
legend(' Weak error',' MC error',2)

if pass==1
  axis([.005 0.5 0.01 50])
  print('-deps2c','barrier_2.eps')
else
  axis([.005 0.5 0.001 5])
  print('-deps2c','barrier_4.eps')
end

figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 0.8]; set(gcf,'pos',pos);
loglog(hh,abs(err3),'b-*',hh,err4,'r-*')
%title('Lookback -- comparison to 2h approximation')
xlabel('h'); ylabel('Error');
legend(' Weak error',' MC error',2)

if pass==1
  axis([.005 0.5 0.02 50])
  print('-deps2c','lookback_2.eps')
else
  axis([.005 0.5 0.002 20])
  print('-deps2c','lookback_4.eps')
end

end

end
