%
% These are the MLMC tests for the paper with Christoph Reisinger
% on the Zakai equation arising from the modelling of credit derivatives
% by Ben Hambly and others
%
% M.B. Giles, C. Reisinger. 'Stochastic finite differences and
% multilevel Monte Carlo for a class of SPDEs in finance',
% SIAM Journal of Financial Mathematics, 3(1):572-592, 2012. 
%

function spde

close all; clear all;

nvert = 3;

N   = 1000;    % samples for convergence tests
L   = 6;       % levels for convergence tests 

N0  = 50;      % initial samples on coarsest levels
Lmin = 2;      % minimum refinement level
Lmax = 6;      % maximum refinement level

Eps = [ 0.001 0.002 0.005 0.01 0.02 ];

fp = fopen('spde.txt','w');
mlmc_test(@spde_l, N,L, N0,Eps,Lmin,Lmax, fp);
fclose(fp);

nvert = 3;      % plotting option (1 for slides, 2 for papers, 3 for full set)

mlmc_plot('spde', nvert);

if(nvert==1)
  figure(1)
  print('-deps2','spde_a.eps')
  figure(2)
  print('-deps2','spde_b.eps')
else
  print('-deps2','spde.eps')
end

%
% Code to solve Christoph's SPDE using a multilevel approach
%

function [sums cost] = spde_l(l,M)

%
% set parameters
%

% old
%
% rho  = 0.5;
% mu   = 0.05;
% x0   = 3.5;

% new

sigma = 0.22;
r     = 0.042;
mu    = (r-0.5*sigma^2)/sigma;
x0    = 5;
rho   = 0.2;

T     = 5;
sqrho = sqrt(rho);

xmax = 16;

%
% calculate number of timesteps and grid intervals
%

% old code based on stability limit
%
% J0 = 4;
% N0 = (2+2*rho)*T*(J0/xmax)^2;
% N0 = ceil(N0);

% new code so N is a multiple of 4*T to give
% a value each quarter

J0 = 10;
N0 = 4*T;

J  = J0 * 2^l;
N  = N0 * 4^l;

xf  = linspace(0,xmax,J+1)';
dxf = xmax / J;
dtf = T / N;
jf  = 2:J;

xc  = linspace(0,xmax,J/2+1)';
dxc = 2*dxf;
dtc = 4*dtf;
jc  = 2:J/2;

sums(1:6) = 0;

%
% main loop, doing M2 at a time
%

M2 = 10;     % max number of samples at a time

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

  Lf = zeros(4*T,m2);
  Lc = zeros(4*T,m2);

  vf = [0; max(0,1-abs(xf(jf)-x0)/dxf)/dxf; 0];      % FE projection of Dirac input at x0
  vc = [0; 0.25*vf(2*jc-2) + 0.5*vf(2*jc-1) + 0.25*vf(2*jc); 0];

  vf = repmat(vf,1,m2);
  vc = repmat(vc,1,m2);

  dMc = 0;

  for n = 1:N
    dMf = sqrt(dtf)*randn(1,m2);
    dMc = dMc + dMf;

    delMf   = repmat(dMf,J-1,1);
    dvdxf   = vf(jf+1,:) - vf(jf-1,:);
    d2vdx2f = vf(jf+1,:) - 2*vf(jf,:) + vf(jf-1,:);

    vf(jf,:) = vf(jf,:) + (0.5*dtf/dxf^2)*((rho/dtf)*delMf.^2+(1-rho)).*d2vdx2f ...
                        - (0.5*dtf/dxf)*(mu + (sqrho/dtf)*delMf).*dvdxf;

    if (mod(n,4)==0)
      delMc   = repmat(dMc,J/2-1,1);
      dvdxc   = vc(jc+1,:) - vc(jc-1,:);
      d2vdx2c = vc(jc+1,:) - 2*vc(jc,:) + vc(jc-1,:);

      vc(jc,:) = vc(jc,:) + (0.5*dtc/dxc^2)*((rho/dtc)*delMc.^2+(1-rho)).*d2vdx2c ...
                          - (0.5*dtc/dxc)*(mu + (sqrho/dtc)*delMc).*dvdxc;

      dMc = 0;
    end

%
% compute fine and coarse path losses at each quarter
%
    R = 0.4;
    if (mod(4*T*n,N)==0)
      Lf((4*T*n)/N,:) = (1-R)*(1 - sum(vf,1)*dxf);
      Lc((4*T*n)/N,:) = (1-R)*(1 - sum(vc,1)*dxc);
    end
  end

%
% calculate discounted tranche loss
%

  a  = 0;
  d  = 0.03;

  Xf = [(d-a)*ones(1,m2); max(d-Lf,0) - max(a-Lf,0)];
  Xc = [(d-a)*ones(1,m2); max(d-Lc,0) - max(a-Lc,0)];

  r    = 0.042;
  n    = 1:4*T;
  disc = diag(exp(-r*n/4));

  Pf = sum(disc*(Xf(n,:)-Xf(n+1,:)),1) / (d-a);
  Pc = sum(disc*(Xc(n,:)-Xc(n+1,:)),1) / (d-a);

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

cost = M*J*N;
% sums(3)/M
