
function adjoint_test

%
% European call
%

S0  = 100;
r   = 0.05;
sig = 0.2;
K   = 60;

T  = 1;
N  = 100;
dt = T/N;

S(1) = S0;
S_1  = 1;   % for delta
S_2  = 0;   % for vega

dW  = sqrt(dt)*randn(1,N);

for n=1:N
  d      = (1 + r*dt + sig*dW(n));
  S_1    = d*S_1;
  S_2    = d*S_2 + S(n)*dW(n);
  S(n+1) = d*S(n);
end

if S(N+1)>K
  P_1 = exp(-r*T) *  S_1;
  P_2 = exp(-r*T) *  S_2;
  P   = exp(-r*T) * (S(N+1)-K);
  Sb  = exp(-r*T);
else
  P_1 = 0;
  P_2 = 0;
  P   = 0;
  Sb  = 0;
end

S0b  = 0;
sigb = 0;

for n=N:-1:1
  sigb = sigb + Sb*S(n)*dW(n);
  Sb   = (1 + r*dt + sig*dW(n))*Sb;
end

S0b = S0b + Sb;

fprintf('\n European call \n');
fprintf(' forward mode sensitivities: delta = %f, vega = %f \n',P_1,P_2);
fprintf(' reverse mode sensitivities: delta = %f, vega = %f \n',S0b,sigb);

%
% re-do using same Brownian increments for CVT check
%

epsi = 1e-10;

for pass=1:2
  if pass == 1
    S(1) = S0 + i*epsi;
  else
    S(1) = S0;
    sig  = sig + i*epsi;
  end

  for n=1:N
    d      = (1 + r*dt + sig*dW(n));
    S(n+1) = d*S(n);
  end

  if real(S(N+1))>K
    P = exp(-r*T) * (S(N+1)-K);
  else
    P = 0;
  end

  if pass == 1
    P_1 = imag(P)/epsi;
  else
    P_2 = imag(P)/epsi;
  end
end

fprintf(' CVT sensitivities:          delta = %f, vega = %f \n',P_1,P_2);

%
% Asian call
%

S0  = 100;
r   = 0.05;
sig = 0.2;
K   = 60;

T  = 1;
N  = 100;
dt = T/N;

S    = [ S0; 0 ];
S_1  = [  1; 0 ];   % for delta
S_2  = [  0; 0 ];   % for vega

dW  = sqrt(dt)*randn(1,N);

for n=1:N
  d      = 1+r*dt+sig*dW(n);
  D      = [ d 0; dt 1 ];
  S_1    = D*S_1;
  S_2    = D*S_2 + [ S(1,n)*dW(n); 0];
  S(:,n+1) = [ d*S(1,n); S(2,n)+dt*S(1,n) ];
end

if S(2,N+1)>K
  P_1 = exp(-r*T) *  S_1(2);
  P_2 = exp(-r*T) *  S_2(2);
  P   = exp(-r*T) * (S(2,N+1)-K);
  Sb  = exp(-r*T) * [ 0; 1 ];
else
  P_1 =  0;
  P_2 =  0;
  P   =  0;
  Sb  = [0; 0];
end

S0b  = 0;
sigb = 0;

for n=N:-1:1
  d    = 1+r*dt+sig*dW(n);
  D    = [ d 0; dt 1 ];
  sigb = sigb + Sb(1)*S(1,n)*dW(n);
  Sb   = D'*Sb;
end

S0b = S0b + Sb(1);

fprintf('\n Asian call \n');
fprintf(' forward mode sensitivities: delta = %f, vega = %f \n',P_1,P_2);
fprintf(' reverse mode sensitivities: delta = %f, vega = %f \n',S0b,sigb);

%
% re-do using same Brownian increments for CVT check
%

for pass=1:2
  if pass==1
    S   = [ S0+i*epsi; 0 ];
  else
    S   = [ S0; 0 ];
    sig = sig + i*epsi;
  end

  for n=1:N
    d      = 1+r*dt+sig*dW(n);
    D      = [ d 0; dt 1 ];
    S(:,n+1) = [ d*S(1,n); S(2,n)+dt*S(1,n) ];
  end

  if real(S(2,N+1))>K
    P = exp(-r*T) * (S(2,N+1)-K);
  else
    P =  0;
  end

  if pass == 1
    P_1 = imag(P)/epsi;
  else
    P_2 = imag(P)/epsi;
  end
end

fprintf(' CVT sensitivities:          delta = %f, vega = %f \n',P_1,P_2);


%
% down-and-out call
%

rng('default');

S0  = 100;
r   = 0.05;
sig = 0.2;
K   = 60;
B   = 90;

T  = 1;
N  = 10;
dt = T/N;

S    = [ S0; 1 ];
S_1  = [  1; 0 ];   % for delta
S_2  = [  0; 0 ];   % for vega

dW  = sqrt(dt)*randn(1,N);

for n=1:N
  d  = 1+r*dt+sig*dW(n);
  p  = exp(- max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / (0.5*sig^2*S(1,n)^2*dt));
  fac1 = (   one(  S(1,n)-B) * max(0,d*S(1,n)-B) ...
         +   max(0,S(1,n)-B) * one(  d*S(1,n)-B) * d ...
         - 2*max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / S(1,n) ) / (0.5*sig^2*S(1,n)^2*dt);
  fac2 = (   max(0,S(1,n)-B) * one(  d*S(1,n)-B) * S(1,n)*dW(n) ...
         - 2*max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / sig    ) / (0.5*sig^2*S(1,n)^2*dt);

  D      = [ d 0; S(2,n)*p*fac1 (1-p) ];
  S_1    = D*S_1;
  S_2    = D*S_2 + [ S(1,n)*dW(n); S(2,n)*p*fac2];
%  fprintf('%f  %f \n',S_1(1),S_1(2))
  S(:,n+1) = [ d*S(1,n); S(2,n)*(1-p) ];
end

if S(1,N+1)>K
  P_1 = exp(-r*T) * ( S_1(1)*S(2,N+1) + (S(1,N+1)-K)*S_1(2) );
  P_2 = exp(-r*T) * ( S_2(1)*S(2,N+1) + (S(1,N+1)-K)*S_2(2) );
  P   = exp(-r*T) * (S(1,N+1)-K)*S(2,N+1);
  Sb  = exp(-r*T) * [ S(2,N+1); (S(1,N+1)-K) ];
else
  P_1 =  0;
  P_2 =  0;
  P   =  0;
  Sb  = [0; 0];
end

S0b  = 0;
sigb = 0;

for n=N:-1:1
  d  = 1+r*dt+sig*dW(n);
  p  = exp(- max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / (0.5*sig^2*S(1,n)^2*dt));
  fac1 = (   one(  S(1,n)-B) * max(0,d*S(1,n)-B) ...
         +   max(0,S(1,n)-B) * one(  d*S(1,n)-B) * d ...
         - 2*max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / S(1,n) ) / (0.5*sig^2*S(1,n)^2*dt);
  fac2 = (   max(0,S(1,n)-B) * one(  d*S(1,n)-B) * S(1,n)*dW(n) ...
         - 2*max(0,S(1,n)-B) * max(0,d*S(1,n)-B) / sig    ) / (0.5*sig^2*S(1,n)^2*dt);

  D      = [ d 0; S(2,n)*p*fac1 (1-p) ];

  sigb = sigb + Sb(1)*S(1,n)*dW(n) + Sb(2)*S(2,n)*p*fac2;
  Sb   = D'*Sb;
end

S0b = S0b + Sb(1);

fprintf('\n down-and-out call \n');
fprintf(' forward mode sensitivities: delta = %f, vega = %f \n',P_1,P_2);
fprintf(' reverse mode sensitivities: delta = %f, vega = %f \n',S0b,sigb);

%
% re-do using same Brownian increments for CVT check
%

for pass=1:2
  if pass==1
    S   = [ S0+i*epsi; 1 ];
  else
    S   = [ S0; 1 ];
    sig = sig + i*epsi;
  end

  for n=1:N
    d = 1+r*dt+sig*dW(n);
    p = exp(- exmax(0,S(1,n)-B) * exmax(0,d*S(1,n)-B) / (0.5*sig^2*S(1,n)^2*dt));
    S(:,n+1) = [ d*S(1,n); S(2,n)*(1-p) ];

    if pass==1
%      fprintf('%f  %f  \n',imag(S(1,n+1))/epsi,imag(S(2,n+1))/epsi)
    end
  end

  if real(S(1,N+1)) > K
    P = exp(-r*T) * (S(1,N+1)-K)*S(2,N+1);
  else
    P =  0;
  end

  if pass == 1
    P_1 = imag(P)/epsi;
  else
    P_2 = imag(P)/epsi;
  end
end

fprintf(' CVT sensitivities:          delta = %f, vega = %f \n',P_1,P_2);

end


function y = one(x)

y = 0.5*(sign(x)+1.0);

end


function z = exmax(x,y)

if real(y) > real(x)
  z = y;
else
  z = x;
end

end

