%
% This tests the use of multilevel Monte Carlo for
% basket options based on 5 underlying assets each
% modelled using Geometric Brownian Motion with
% correlation between the driving Brownian motions
%

function basket

close all; clear all;

global option

nvert = 1;
M     = 2;      % refinement cost factor
 
for option = 1:5

  if (option==1) 
    stitle = 'GBM model: European call';
    N      = 100000;   % samples for convergence tests
    L      = 8;        % levels for convergence tests 
    N0     = 100;      % initial samples on first 3 levels
    Eps    = [ 0.01 0.02 0.05 0.1 0.2 ];

    figs = mlmc_test(@basket_l, M, N,L, N0,Eps, nvert);

    subplot(nvert,2,2*nvert-1); axis([0 8 50 1e7])
    legend('0.01','0.02','0.05','0.1','0.2',1)
    subplot(nvert,2,2*nvert); axis([1e-2 2e-1 5e2 5e5])
    figure(figs(1))
    subplot(nvert,2,1); axis([0 8 -17.5 12.5])
    subplot(nvert,2,2); axis([0 8 -12.5 7.5])

  elseif (option==2) 
    stitle = 'GBM model: Asian option';
    N      = 100000;   % samples for convergence tests
    L      = 8;        % levels for convergence tests 
    N0     = 100;      % initial samples on each level
    Eps    = [ 0.01 0.02 0.05 0.1 0.2 ];

    figs = mlmc_test(@basket_l, M, N,L, N0,Eps, nvert);

    subplot(nvert,2,2*nvert-1); axis([0 8 1e2 1e8])
    legend('0.01','0.02','0.05','0.1','0.2',1)
    subplot(nvert,2,2*nvert); axis([1e-2 2e-1 2e2 5e5])
    figure(figs(1))
    subplot(nvert,2,1); axis([0 8 -17.5 12.5])
    subplot(nvert,2,2); axis([0 8 -12.5 7.5])

  elseif (option==3)
    stitle = 'GBM model: lookback option';
    N      = 100000;   % samples for convergence tests
    L      = 8;        % levels for convergence tests 
    N0     = 100;      % initial samples on each level
    Eps    = [ 0.01 0.02 0.05 0.1 0.2 ];

    figs = mlmc_test(@basket_l, M, N,L, N0,Eps, nvert);

    subplot(nvert,2,2*nvert-1); axis([0 8 1e2 1e8])
    legend('0.01','0.02','0.05','0.1','0.2',1)
    subplot(nvert,2,2*nvert); axis([1e-2 2e-1 5e2 2e6])
    figure(figs(1))
    subplot(nvert,2,1); axis([0 8 -17.5 12.5])
    subplot(nvert,2,2); axis([0 8 -12.5 7.5])

  elseif (option==4) 
    stitle = 'GBM model: digital call';
    N      = 500000;   % samples for convergence tests
    L      = 8;        % levels for convergence tests 
    N0     = 200;      % initial samples on each level
    Eps    = [ 0.01 0.02 0.05 0.1 0.2 ];

    figs = mlmc_test(@basket_l, M, N,L, N0,Eps, nvert);

    subplot(nvert,2,2*nvert-1); axis([0 8 1e2 1e8])
    legend('0.01','0.02','0.05','0.1','0.2',1)
    subplot(nvert,2,2*nvert); axis([1e-2 2e-1 2e2 1e7])
    figure(figs(1))
    subplot(nvert,2,1); axis([0 8 -17.5 12.5])
    subplot(nvert,2,2); axis([0 8 -12.5 7.5])

  elseif (option==5) 
    stitle = 'GBM model: barrier call';
    N      = 500000;   % samples for convergence tests
    L      = 8;        % levels for convergence tests 
    N0     = 100;      % initial samples on each level
    Eps    = [ 0.01 0.02 0.05 0.1 0.2 ];

    figs = mlmc_test(@basket_l, M, N,L, N0,Eps, nvert);

    subplot(nvert,2,2*nvert-1); axis([0 8 1e2 1e8])
    legend('0.01','0.02','0.05','0.1','0.2',1)
    subplot(nvert,2,2*nvert); axis([1e-2 2e-1 5e2 2e6])
    figure(figs(1))
    subplot(nvert,2,1); axis([0 8 -17.5 12.5])
    subplot(nvert,2,2); axis([0 8 -12.5 7.5])
  end

%
% print figures
%
  if(nvert==1)
    figure(figs(1))
    print('-deps2',sprintf('basket%da.eps',option))
    figure(figs(2))
    print('-deps2',sprintf('basket%db.eps',option))
  else
    print('-deps2',sprintf('basket%d.eps',option))
  end
end

%-------------------------------------------------------
%
% level l estimator
%

function [sum1 sum2] = basket_l(l,N)

global option

dim = 5;
Sig = eye(dim) +  0.25*(ones(dim)-eye(dim));
C   = chol(Sig, 'lower');
alf = ones(dim,1)/dim;


K   = 100;

T   = 1;
r   = 0.05;
sig = diag(0.2+0.05*(1:dim));

B   = 0.85*K;

nf = 2^l;
nc = nf/2;


hf = T/nf;
hc = T/nc;

if option==4
  nc = nc-1;
end

sum1(1:4) = 0;
sum2(1:2) = 0;

for N1 = 1:10000:N
  N2 = min(10000,N-N1+1);

%
% multivariate GBM model
%

  alfm = repmat(alf,1,N2);

  X0 = K;  % start at strike

  Xf = X0*ones(dim,N2);
  Xc = Xf;

  Af  = 0.5*hf*alf'*Xf;
  Ac  = 0.5*hc*alf'*Xc;

  Mf  = alf'*Xf;
  Mc  = alf'*Xc;

  Bf  = 1;
  Bc  = 1;

  if l==0 & option~=4
    dWf = sqrt(hf)*C*randn(dim,N2);
    dIf = sqrt(hf/12)*hf*C*randn(dim,N2);

    Xf0 = Xf;
    Xf  = Xf + r*Xf*hf + sig*Xf.*dWf + 0.5*sig^2*(Xf.*(dWf.^2-hf));

    vf  = sig*Xf0;
    Af  = Af + 0.5*hf*alf'*Xf + alf'*(vf.*dIf);

    Lf  = log(rand(1,N2));

    vf2 = sum( (alfm.*vf).* (Sig*(alfm.*vf)) );
    Mf  = min(Mf,0.5*(alf'*(Xf0+Xf)-sqrt((alf'*(Xf-Xf0)).^2-2*hf*vf2.*Lf)));
    Bf  = Bf.*(1-exp(-2*max(0,(alf'*Xf0-B).*(alf'*Xf-B)./(hf*vf2))));

  else
    for n = 1:nc
      Lf = log(rand(2,N2));

      for m = 1:2
        dWf = sqrt(hf)*C*randn(dim,N2);
        dIf = sqrt(hf/12)*hf*C*randn(dim,N2);

        Xf0 = Xf;
        Xf  = Xf + r*Xf*hf + sig*Xf.*dWf + 0.5*sig^2*(Xf.*(dWf.^2-hf));

        vf  = sig*Xf0;
        Af  = Af + hf*alf'*Xf + alf'*(vf.*dIf);

        vf2 = sum( (alfm.*vf).* (Sig*(alfm.*vf)) );
        Mf  = min(Mf,0.5*(alf'*(Xf0+Xf)-sqrt((alf'*(Xf-Xf0)).^2-2*hf*vf2.*Lf(m,:))));
        Bf  = Bf.*(1-exp(-2*max(0,(alf'*Xf0-B).*(alf'*Xf-B)./(hf*vf2))));

        if m==1
          dWc = dWf;
          ddW = dWf;
          dIc = dIf;
        else
          dWc = dWc + dWf;
          ddW = ddW - dWf;
          dIc = dIc + dIf;
        end
      end

      Xc0 = Xc;
      Xc  = Xc + r*Xc*hc + sig*Xc.*dWc + 0.5*sig^2*Xc.*(dWc.^2-hc);

      vc  = sig*Xc0;
      Ac  = Ac + hc*alf'*Xc + alf'*(vc.*(dIc+0.25*hc*ddW));

      vc2 = sum( (alfm.*vc).* (Sig*(alfm.*vc)) );

      Xc1 = 0.5*(Xc0 + Xc + vc.*ddW);
      Mc  = min(Mc, 0.5*(alf'*(Xc0+Xc1)-sqrt((alf'*(Xc1-Xc0)).^2-2*hf*vc2.*Lf(1,:))));
      Mc  = min(Mc, 0.5*(alf'*(Xc1+Xc) -sqrt((alf'*(Xc -Xc1)).^2-2*hf*vc2.*Lf(2,:))));

      Bc  = Bc .*(1-exp(-2*max(0,(alf'*Xc0-B).*(alf'*Xc1-B)./(hf*vc2))));
      Bc  = Bc .*(1-exp(-2*max(0,(alf'*Xc1-B).*(alf'*Xc -B)./(hf*vc2))));
    end

    Af = Af - 0.5*hf*alf'*Xf;
    Ac = Ac - 0.5*hc*alf'*Xc;
  end

  if option==1
    Pf  = max(0,alf'*Xf-K);
    Pc  = max(0,alf'*Xc-K);
  elseif option==2
    Pf  = max(0,Af-K);
    Pc  = max(0,Ac-K);
  elseif option==3
    Pf  = alf'*Xf - Mf;
    Pc  = alf'*Xc - Mc;
  elseif option==4
    if(l==0)
      vf  = sig*Xf;
      vf2 = sum( (alfm.*vf).* (Sig*(alfm.*vf)) );
      Pf  = 100*ncf((alf'*(Xf+r*Xf*hf)-K)./sqrt(hf*vf2));
      Pc  = 100*Pf;
    else
      vf  = sig*Xf;
      vf2 = sum( (alfm.*vf).* (Sig*(alfm.*vf)) );
      vc  = sig*Xc;
      vc2 = sum( (alfm.*vc).* (Sig*(alfm.*vc)) );
      dWf = sqrt(hf)*C*randn(dim,N2);
      Xf  = Xf + r*Xf*hf + vf.*dWf + 0.5*sig^2*Xf.*(dWf.^2-hf);
      vf  = sig*Xf;
      vf2 = sum( (alfm.*vf).* (Sig*(alfm.*vf)) );
      Pf  = 100*ncf((alf'*(Xf+r*Xf*hf)-K)./sqrt(hf*vf2));
      Pc  = 100*ncf((alf'*(Xc+r*Xc*hc+vc.*dWf)-K)./sqrt(hf*vc2));
    end
  elseif option==5
    Pf  = Bf.*max(0,alf'*Xf-K);
    Pc  = Bc.*max(0,alf'*Xc-K);
  end

  dP  = exp(-r*T)*(Pf-Pc);
  Pf  = exp(-r*T)*Pf;

  if l==0
    dP = Pf;
  end

  sum1(1) = sum1(1) + sum(dP);
  sum1(2) = sum1(2) + sum(dP.^2);
  sum1(3) = sum1(3) + sum(dP.^3);
  sum1(4) = sum1(4) + sum(dP.^4);
  sum2(1) = sum2(1) + sum(Pf);
  sum2(2) = sum2(2) + sum(Pf.^2);
end

