%
% These are similar to the MLMC tests for basket options
% shown at the Global Derivatives conference and in the
% 2009 Winter Simulation Conference paper
%
% However, these use MLQMC
%

function basket

close all; clear all;

global option

addpath('..');

M     = 2;      % refinement cost factor

Nr1   = 512;    % number of randomisations for initial convergence tests
%N     = [1 16 256 4096]; % QMC set sizes for initial convergence tests
N     = [1 16 256]; % QMC set sizes for initial convergence tests

Nr    = 32;     % number of randomisations for main MLQMC calcs
Lmin  = 2;      % minimum refinement level
Lmax  = 10;     % maximum refinement level

Qdims = 5 * M.^(0:Lmax); % QMC dimension on each level

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

  filename = ['basket_' num2str(option)];
  fp = fopen([filename '.txt'],'w');
  mlqmc_test(@basket_l, M,Qdims, N,Nr1,L, Nr,Eps,Lmin,Lmax, fp);
  fclose(fp);

%
% plot results
%
  nvert = 3;
  mlqmc_plot(filename, nvert);

  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 [dP Pf] = basket_l(l,dimUQ,N,UQ)

global option

dim = 5;

if dimUQ ~= dim*2^l
  error('wrong QMC dimension')
end

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

%
% multivariate GBM model
%

alfm = repmat(alf,1,N);

X0 = K;  % start at strike

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

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

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

Bf  = 1;
Bc  = 1;

Z  = ncfinv(UQ);    % convert quasi-uniforms into quasi-Normals
dW = bb(Z,T,C);     % convert quasi-Normals into correlated Brownian
                    % increments using Brownian Bridge construction

if l==0 & option~=4
  Lf  = log(rand(1,N));
  dWf = dW;
  dIf = sqrt(hf/12)*hf*C*randn(dim,N);

  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);

  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,N));

    for m = 1:2
      dWf = dW((2*n+m-3)*dim+1:(2*n+m-2)*dim,:);
      dIf = sqrt(hf/12)*hf*C*randn(dim,N);

      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));
  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,N);
    dWf = dW((2*nc)*dim+1:(2*nc+1)*dim,:);
    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

if (l==0)
  Pc = zeros(1,N);
end

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

