%
% function V = heston_call(r,kappa,theta,omega,rho,V0,T,S0,K,opt,opt2)
%
% analytic solution for European call option in Heston model, using
% algorithm in "Not-so-complex logarithms in the Heston model"
% by Christian Kahl and Peter Jaekel
% http://www.math.uni-wuppertal.de/~kahl/publications/NotSoComplexLogarithmsInTheHestonModel.pdf
% http://www.btinternet.com/~pjaeckel/NotSoComplexLogarithmsInTheHestonModel.pdf
%
% r     - interest rate
% kappa - mean reversion rate
% theta - mean reversion volatility
% omega - vol-of-vol coefficient
% rho   - correlation factor
% V0    - initial volatility
% T     - time interval
% S     - initial asset price
% K     - strike price
% opt   - 'value', 'delta' or 'gamma'
% opt2  - (optional) '1' (default) or '2'
%         '1' uses complex arithmetic (unsuitable for CVT check)
%         '2' uses real arithmetic (with complex extension for CVT check)
% V     - option value

function V = heston_call(r,kappa,theta,omega,rho,V0,T,S,K,opt,opt2)

if nargin == 10
  opt2 = '1';
elseif nargin ~= 11
  error('wrong number of arguments');
end

if opt2=='1'
  V = quadl(@(x)integrand1(x,r,kappa,theta,omega,rho,V0,T,S,K,opt),0,1,1e-12);
else
  V = quadl(@(x)integrand2(x,r,kappa,theta,omega,rho,V0,T,S,K,opt),0,1,1e-12);
end

V = exp(-r*T)*V;

%--------------------------------------------------------

function y = integrand1(x,r,kappa,theta,omega,rho,V0,T,S,K,opt)

x    = max(1e-20,min(x,1-1e-10));
Cinf = sqrt(1-rho^2)*(V0+kappa*theta*T)/omega;
u    = - log(x)/Cinf;
F    = S*exp(r*T);

for ipass = 1:2
  um  = u - i*(2-ipass);
  d   = sqrt( (rho*omega*um*i - kappa).^2 + omega^2*(um*i + um.^2));
  c   = (kappa-rho*omega*um*i+d)./(kappa-rho*omega*um*i-d);
  tc  = angle(c);
  GD  = c-1;
  m   = floor( (tc+pi)/(2*pi) );
  GN  = c.*exp(i*imag(d)*T) - exp(-real(d)*T);
  n   = floor( (tc+imag(d)*T+pi)/(2*pi) );
  lnG = real(d)*T + log(abs(GN)./abs(GD)) + i*(angle(GN)-angle(GD)+2*pi*(n-m));
  D   = ((kappa-rho*omega*um*i+d)/omega^2) .* ((exp(d*T)-1)./(c.*exp(d*T)-1));
  C   = ((kappa*theta)/omega^2) * ( (kappa-rho*omega*um*i+d)*T - 2*lnG );
  phi = exp(C+D*V0+i*um*log(F));

  switch lower(opt)
    case 'value'
      f(ipass,:) = real( exp(-i*u*log(K)).*phi ./ (i*u) );

    case 'delta'
      f(ipass,:) = real( um.*exp(-i*u*log(K)).*phi ./ (u*S) );

    case 'gamma'
      f(ipass,:) = real( um.*(i*um-1).*exp(-i*u*log(K)).*phi ./ (u*S^2) );

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

y = (f(1,:) - K*f(2,:))./(x*pi*Cinf);

switch lower(opt)
  case 'value'
    y = y + 0.5*(F-K);

  case 'delta'
    y = y + 0.5*exp(r*T);
end

%x
%y


%--------------------------------------------------------

function y = integrand2(x,r,kappa,theta,omega,rho,V0,T,S,K,opt)

x    = max(1e-20,min(x,1-1e-10));
Cinf = sqrt(1-rho^2)*(V0+kappa*theta*T)/omega;
u    = - log(x)/Cinf;
F    = S*exp(r*T);

for ipass = 1:2
  um_r = u;
  um_i = - (2-ipass);

  s_r  = -rho*omega*um_i - kappa;
  s_i  =  rho*omega*um_r;

  [t1_r, t1_i] = cprod( s_r, s_i,  s_r, s_i);
  [t2_r, t2_i] = cprod(um_r,um_i, um_r,um_i);

  d2_r = t1_r + omega^2*(t2_r - um_i);
  d2_i = t1_i + omega^2*(t2_i + um_r);

  [d_r, d_i] = csqrt(d2_r, d2_i);
  [c_r, c_i] = cdiv(kappa+rho*omega*um_i+d_r,-rho*omega*um_r+d_i, ...
                    kappa+rho*omega*um_i-d_r,-rho*omega*um_r-d_i);
  tc = catan2(c_i,c_r);

  GD_r = c_r-1;
  GD_i = c_i;

  m    = floor( (real(tc)+pi)/(2*pi) );

  GN_r = c_r.*cos(d_i*T) - c_i.*sin(d_i*T) - exp(-d_r*T);
  GN_i = c_r.*sin(d_i*T) + c_i.*cos(d_i*T);

  n    = floor( (real(tc)+d_i*T+pi)/(2*pi) );

  lnG_r = d_r*T + 0.5*log((GN_r.^2+GN_i.^2)./(GD_r.^2+GD_i.^2));
  lnG_i = catan2(GN_i,GN_r)-catan2(GD_i,GD_r)+2*pi*(n-m);

  [p0_r, p0_i] = cexp(d_r*T, d_i*T);
  [p1_r, p1_i] = cprod(c_r,c_i, p0_r,p0_i);
  [p2_r, p2_i] = cdiv(p0_r-1,p0_i, p1_r-1, p1_i);
  [D_r, D_i]   = cprod(kappa+rho*omega*um_i+d_r, -rho*omega*um_r+d_i, ...
                       p2_r/omega^2, p2_i/omega^2);

  C_r = ((kappa*theta)/omega^2) * ( (kappa+rho*omega*um_i+d_r)*T - 2*lnG_r );
  C_i = ((kappa*theta)/omega^2) * ( (     -rho*omega*um_r+d_i)*T - 2*lnG_i );

  [phi_r,phi_i] = cexp(C_r+D_r*V0-um_i*log(F), C_i+D_i*V0+um_r*log(F));

  [p0_r, p0_i] = cprod(cos(u*log(K)),-sin(u*log(K)), phi_r,phi_i);

  switch lower(opt)
    case 'value'
      f(ipass,:) = p0_i./u;

    case 'delta'
      [p1_r, p1_i] = cprod(p0_r,p0_i, um_r,um_i);
      f(ipass,:)   = p1_r./(u.*S);

    case 'gamma'
      [p1_r, p1_i] = cprod(p0_r,p0_i, um_r,um_i);
      [p1_r, p1_i] = cprod(p1_r,p1_i, um_r,um_i+1);
      f(ipass,:)   = - p1_i./(u.*S.^2);

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

y = (f(1,:) - K*f(2,:))./(x*pi*Cinf);

switch lower(opt)
  case 'value'
    y = y + 0.5*(F-K);

  case 'delta'
    y = y + 0.5*exp(r*T);
end

%x
%y


%
% complex functions implemented using real arithmetic
%

function [c_r, c_i] = cprod(a_r,a_i,b_r,b_i)

c_r = a_r.*b_r - a_i.*b_i;
c_i = a_r.*b_i + a_i.*b_r;


function [c_r, c_i] = cdiv(a_r,a_i,b_r,b_i)

d_r =   b_r ./ (b_r.^2+b_i.^2);
d_i = - b_i ./ (b_r.^2+b_i.^2);

c_r = a_r.*d_r - a_i.*d_i;
c_i = a_r.*d_i + a_i.*d_r;


function [c_r, c_i] = csqrt(a_r,a_i)

c_r = sqrt(0.5*(a_r + sqrt(a_r.^2+ a_i.^2)));
c_i = 0.5*a_i./c_r;


function [c_r, c_i] = cexp(a_r,a_i)

c_r = exp(a_r).*cos(a_i);
c_i = exp(a_r).*sin(a_i);

%
% complex extension to atan2 function
%

function c = catan2(a,b)

c = atan2(real(a),real(b)) ...
  + (a./b - real(a)./real(b))./(1+a.^2./b.^2);


%
% first check consistency between integrand1/2
%

function integrand_check

r     = 0.05;
kappa = 5;
theta = 0.04;
omega = 0.25;
rho   = -0.5;
S     = 1;
K     = 1.2;
V0    = 0.04;
T     = 0.5;

for iopt = 1:3
  if iopt ==1
    opt = 'value'
  elseif iopt ==2
    opt = 'delta'
  elseif iopt ==3
    opt = 'gamma'
   end
  integrand1(0.1,r,kappa,theta,omega,rho,V0,T,S,K,opt)
  integrand2(0.1,r,kappa,theta,omega,rho,V0,T,S,K,opt)
end

%
% second check correctness of integrand2 using CVT 
%

eps = 1e-10;

delta = integrand2(0.1,r,kappa,theta,omega,rho,V0,T,S,K,'delta');
V     = integrand2(0.1,r,kappa,theta,omega,rho,V0,T,S+i*eps,K,'value');

disp(sprintf('delta (analytic) = %g',delta))
disp(sprintf('delta (CVT)      = %g',imag(V)/eps))

gamma = integrand2(0.1,r,kappa,theta,omega,rho,V0,T,S,K,'gamma');
delta = integrand2(0.1,r,kappa,theta,omega,rho,V0,T,S+i*eps,K,'delta');

disp(sprintf('gamma (analytic) = %g',gamma))
disp(sprintf('gamma (CVT)      = %g',imag(delta)/eps))

figure
x  = linspace(0,1,10000);
y1 = integrand1(x,r,kappa,theta,omega,rho,V0,T,S,K,'value');
y2 = integrand2(x,r,kappa,theta,omega,rho,V0,T,S,K,'value');
norm(y1-y2)
plot(x,y1,x,y2);
axis([0 .1 -10 10])


