%
% These are the MLMC tests for basket options shown
% at the Global Derivatives conference and in the
% 2009 Winter Simulation Conference paper
%
% There are some minor differences:
% -- the new MLMC driver is a little different
%    (slightly more efficient and more robust)
% -- MATLAB has changed the default RNG
%

function basket

close all; clear all;

addpath('..');

N0    = 200;    % initial samples on coarse levels
Lmin  = 2;      % minimum refinement level
Lmax  = 10;     % maximum refinement level

for option = 1:5
  if (option==1) 
    fprintf(1,'\n ---- European call ---- \n');
    N    = 50000;    % samples for convergence tests
    L    = 8;        % levels for convergence tests 
    Eps  = [ 0.01 0.02 0.05 0.1 0.2 ];
  elseif (option==2) 
    fprintf(1,'\n ---- Asian call ---- \n');
    N    = 50000;    % samples for convergence tests
    L    = 8;        % levels for convergence tests 
    Eps  = [ 0.01 0.02 0.05 0.1 0.2 ];
  elseif (option==3)
    fprintf(1,'\n ---- lookback call ---- \n');
    N    = 50000;    % samples for convergence tests
    L    = 8;        % levels for convergence tests 
    Eps  = [ 0.01 0.02 0.05 0.1 0.2 ];
  elseif (option==4) 
    fprintf(1,'\n ---- digital call ---- \n');
    N    = 200000;   % samples for convergence tests
    L    = 8;        % levels for convergence tests 
    Eps  = [ 0.01 0.02 0.05 0.1 0.2 ];
  elseif (option==5) 
    fprintf(1,'\n ---- barrier call ---- \n');
    N    = 200000;   % samples for convergence tests
    L    = 8;        % levels for convergence tests 
    Eps  = [ 0.01 0.02 0.05 0.1 0.2 ];
  end

  filename = ['basket_' num2str(option)];
  fp = fopen([filename '.txt'],'w');
  mlmc_test(@basket_l, N,L, N0,Eps,Lmin,Lmax, fp, option);
  fclose(fp);

%
% plot results
%
  nvert = 2;
  mlmc_plot(filename, nvert);
  
  subplot(nvert,2,2*nvert-1)
  current = axis; axis([current(1:2) 10 current(4)])

  if option==4
    subplot(nvert,2,1)
    current = axis; axis([current(1:2) -10 15])
    subplot(nvert,2,2*nvert-1)
    current = axis; axis([current(1:2) 100 1e8])
    subplot(nvert,2,2*nvert)
    current = axis; axis([current(1:2) 100 1e7])
  end
    
  if(nvert==1)
    figure(1)
    print('-deps2',[filename 'a.eps'])
    figure(2)
    print('-deps2',[filename 'b.eps'])
  else
    print('-deps2',[filename '.eps'])
  end
end

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

function [sums, cost] = basket_l(l,N, option)

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

asian_switch = 1;  % switch for which Asian treatment to use

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

sums(1:6) = 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
    Lf  = log(rand(1,N2));
    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 + asian_switch*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)));
    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 + asian_switch*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 + asian_switch*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  = K*ncf((alf'*(Xf+r*Xf*hf)-K)./sqrt(hf*vf2));
      Pc  = K*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

  sums(1) = sums(1) + sum(dP);
  sums(2) = sums(2) + sum(dP.^2);
  sums(3) = sums(3) + sum(dP.^3);
  sums(4) = sums(4) + sum(dP.^4);
  sums(5) = sums(5) + sum(Pf);
  sums(6) = sums(6) + sum(Pf.^2);
end

cost = N*nf;   % cost defined as number of fine timesteps

