%
% function V = lookback_call(r,sigma,T,S,opt)
%
% Black-Scholes floating strike lookback call option 
% solution as defined in
% https://en.wikipedia.org/wiki/Lookback_option
% with initial S_min equal to initial S
%
% r     - interest rate
% sigma - volatility
% T     - time interval
% S     - initial asset price
% opt   - 'value', 'delta', 'gamma' or 'vega'
% V     - option value(s)
%

function V = lookback_call(r,sigma,T,S,opt)

if nargin ~= 5
  error('wrong number of arguments');
end

d1 = (r+0.5*sigma^2)*T / (sigma*sqrt(T));
d2 = (r-0.5*sigma^2)*T / (sigma*sqrt(T));

k = sigma^2/(2*r);

switch opt
  case 'value'
    V = S.* ( N(d1) - N(-d1)*k - exp(-r*T)*(N(d2) - N(d2)*k) );

  case 'delta'
    V = N(d1) - N(-d1)*k - exp(-r*T)*(N(d2) - N(d2)*k);

  case 'gamma'
    V = 0;

  case 'vega'
    V =  S.*(         (exp(-0.5*d1.^2)./sqrt(2*pi)).*( sqrt(T)-d1/sigma)*(1+k) ...
          - exp(-r*T)*(exp(-0.5*d2.^2)./sqrt(2*pi)).*(-sqrt(T)-d2/sigma)*(1-k) ...
          - (N(-d1)- exp(-r*T)*N(d2))*sigma/r );

  otherwise
    error('invalid value for opt -- must be ''value'', ''delta'', ''gamma'', ''vega''')
end



%
% Normal cumulative distribution function
%

function ncf = N(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

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