
function discrete_barrier

close all;
clear all;

r   = 0.02;
sig = 0.2;
T   = 1;
S0  = 100;
B   = 90;
K   = 100;

N   = 1e6;

for pass = 1:2

  vals  = [];
  vars  = [];
  bumps = [];

  for bump = 10.^linspace(-3,1,20)
    if (pass==1)
      fm = EM_1(r,sig,T,S0-bump,B,K, N);
      fp = EM_1(r,sig,T,S0+bump,B,K, N);
    else
      fm = EM_2(r,sig,T,S0-bump,B,K, N);
      fp = EM_2(r,sig,T,S0+bump,B,K, N);
    end

    Del = (fp-fm)/(2*bump);

    val =  sum(Del)/N;
    var = (sum(Del.^2)/N - val^2)/(N-1);

    vals  = [vals  val];
    vars  = [vars  var];
    bumps = [bumps bump];
  end

  figure; pos=get(gcf,'pos'); pos(3:4)=pos(3:4).*[0.8 1.2]; set(gcf,'pos',pos);
  subplot(2,1,1)
  semilogx(bumps,vals);
  xlabel('bump'); ylabel('Delta')

  subplot(2,1,2)
  loglog(bumps,vars);
  xlabel('bump'); ylabel('variance')

  print('-deps2c',sprintf('results_%1d.eps',pass))
end

%
% complex variable trick for pathwise sensitivity
%

bump = 1e-20;
f   = EM_2(r,sig,T,S0+i*bump,B,K, N);
Del = imag(f)/bump;

val =  sum(Del)/N;
var = (sum(Del.^2)/N - val^2)/(N-1);

fprintf('Delta =    %g \n',val)
fprintf('Variance = %g \n',var)

end



function f = EM_1(r,sig,T,S0,B,K, N)

rng('default');  % reset random number generator

h = 0.05;
t = 0;

S     = S0*ones(1,N);
S_min = S;

for n = 1:20
  t  = t + h;
  dW = sqrt(h)*randn(1,N);
  S  = S.*( 1 + r*h + sig*dW);

  if(mod(n,4)==0)
    S_min = min(S_min,S);
  end
end

f = exp(-r*T)*max(S-K,0) .* (0.5*(sign(S_min-B)+1));

end


function f = EM_2(r,sig,T,S0,B,K, N)

rng('default');  % reset random number generator

h = 0.05;
t = 0;

S     = S0*ones(1,N);
S_min = S;
prob  = ones(1,N);

for n = 1:16
  if ( (n>1) & n<16 & mod(n,3)==1 )  % jump-over
    S_old = S;
    sd    = sqrt(0.5*h)*sig*S;

    t     = t + 2*h;
    dW = sqrt(2*h)*randn(1,N);
    S  = S.*( 1 + 2*r*h + sig*dW);

    S_mid = 0.5*(S_old+S);
    prob = prob .* ncf((S_mid-B)/sd);
  else
    t  = t + h;
    dW = sqrt(h)*randn(1,N);
    S  = S.*( 1 + r*h + sig*dW);
  end
end

f = exp(-r*T)*exmax(S-K,0) .* prob;

end


%
% Normal cumulative distribution function
%

function val = ncf(x)

%ncf = 0.5*(1+erf(x/sqrt(2)));

xr = real(x);
xi = imag(x);

if abs(xi)>1e-10
  error 'imag(x) too large in N(x)'
end

val = 0.5*(1+erf(xr/sqrt(2))) ...
    + i*xi.*exp(-0.5*xr.^2)/sqrt(2*pi);

end

%
% analytic extension of real max function
%

function umax = exmax(u1,u2)

sw   = real(u1) > real(u2);
umax = sw.*u1 + (1-sw).*u2;

end
