function y = HestonFFTVanilla(phi,S,K,T,rd,rf,kappa,theta,sigma,rho,v0,alpha,method)
%HESTONFFTVANILLA European FX option price in the Heston model 
%(Carr-Madan approach).
%   Y = HESTONFFTVANILLA(PHI,S,K,T,R,RF,KAPPA,THETA,SIGMA,RHO,V0) returns
%   the price of a European call (PHI=1) or put (PHI=-1) option given 
%   spot price S, strike K, time to maturity (in years) T, domestic R and 
%   foreign RF interest rates, rate of mean reversion KAPPA, average level 
%   of volatility THETA, volatility of volatility SIGMA, correlation 
%   between the Wiener increments driving the spot and vol processes RHO, 
%   and initial volatility VO.
%   Y = HESTONFFTVANILLA(...,ALPHA,METHOD) allows to specify the coefficient 
%   ALPHA of the exponential smoother (default: ALPHA=0.75 for calls and 
%   1+ALPHA=1.75 for puts) and the integration method of the complex
%   integral: METHOD = 0 (default) -> adaptive Gauss-Kronrod quadrature,
%   or METHOD = 1 -> FFT + Simpson's rule (as in [2]).

clear all;
clc;
format long;

% if (nargin < 11)
%   error ('Wrong number of input arguments.')
% else
%   if (nargin < 13)
%     % Set default value of method: 
%     method = 0; % adaptive Gauss-Kronrod quadrature
%     if (nargin == 11)
%       % Set parameter alpha (see [3])
%       alpha = 0.75; % used for call option
%     end
%   end
% end
% 
% if (phi==-1), %put option
%   alpha = alpha + 1; 
% end

daten=xlsread('CallPreise.xlsx','Jan','B2:K55'); 
phi = 1;
S=daten(:,1);
K=daten(:,2);
T=daten(:,3);
rd=daten(:,4);
rf=daten(:,5);
kappa=daten(:,6);
theta=daten(:,7);
sigma=daten(:,8);
rho=daten(:,9);
v0=daten(:,10);
alpha=0.75;
method=1;

s0 = log(S);
k = log(K);

for count = 1:51

    if (method == 0)
      % Integrate using adaptive Gauss-Kronrod quadrature
      y = exp(-phi.*k(count).*alpha).*quadgk(@(v) HestonFFTVanillaInt(phi,s0(count),k(count),T(count),rd(count),rf(count),kappa(count),theta(count),sigma(count),rho(count),v0(count),alpha,v),0,inf,'RelTol',1e-8)./pi;
        else
          % FFT with Simpson's rule (as suggested in [2])
          N = 2^10; 
          eta = 0.25;
          v =(0:N-1)*eta;

          lambda = 2*pi/(N*eta);
          b = N*lambda/2;
          ku = -b+lambda.*(0:N-1);
          u = v - (phi.*alpha+1)*1i;
          d = sqrt((rho(count).*sigma(count).*u*1i-kappa(count)).^2+sigma(count).^2.*(1i*u+u.^2));
          g = (kappa(count)-rho(count).*sigma(count).*1i.*u-d)./(kappa(count)-rho(count).*sigma(count)*1i.*u+d);

          % Characteristic function (see [1])
          A = 1i*u.*(s0(count) + (rd(count)-rf(count)).*T(count));
          B = theta(count).*kappa(count).*sigma(count).^(-2).*((kappa(count)-rho(count).*sigma(count)*1i.*u-d).*T(count)-2*log((1-g.*exp(-d.*T(count)))./(1-g)));
          C = v0(count).*sigma(count).^(-2).*(kappa(count)-rho(count).*sigma(count)*1i.*u-d).*(1-exp(-d.*T(count)))./(1-g.*exp(-d.*T(count)));
          charFunc = exp(A + B + C); 
          F = charFunc*exp(-rd(count)*T(count))./(alpha^2 + phi.*alpha - v.^2 + 1i*(phi.*2*alpha +1).*v);

          % Use Simpson's approximation to calculate FFT (see [2])
          SimpsonW = 1/3*(3 + (-1).^[1:N] - [1, zeros(1,N-1)]);
          FFTFunc = exp(1i*b*v).*F*eta.*SimpsonW;
          payoff = real(fft(FFTFunc));
          OptionValue = exp(-phi.*ku.*alpha).*payoff./pi;
          % Interpolate to get option price for a given strike
          y = interp1(ku,OptionValue,k);
    end
end

%%%%%%%%%%%%% INTERNALLY USED ROUTINE %%%%%%%%%%%%%


function payoff = HestonFFTVanillaInt(phi,s0,k,T,rd,rf,kappa,theta,sigma,rho,v0,alpha,v)
%HESTONFFTVANILLAINT Auxiliary function used by HESTONFFTVANILLA.
%   PAYOFF=HESTONFFTVANILLAINT(phi,s0,k,T,r,rf,kappa,theta,sigma,rho,v0,alpha,v)
%   returns the values of the auxiliary function evaluated at points V, 
%   given log(spot price) S0, log(strike) K, time to maturity (in years) T,
%   domestic interest rate R, foreign interest rate RF, level of mean
%   reversion KAPPA, long-run variance THETA, vol of vol SIGMA, correlation 
%   RHO, initial volatility VO, damping coefficient ALPHA and option type PHI:
%   PHI = 1 --> call option,
%   PHI = -1 --> put option.
for count = 1:1
   
u = v - (phi.*alpha+1)*1i;
d = sqrt((rho(count).*sigma(count).*u*1i-kappa(count)).^2+sigma(count).^2.*(1i*u+u.^2));
g = (kappa(count)-rho(count).*sigma(count).*1i.*u-d)./(kappa(count)-rho(count).*sigma(count)*1i.*u+d);

% Characteristic function (see [1])
A = 1i*u.*(s0(count) + (rd(count)-rf(count)).*T);
B = theta(count).*kappa(count).*sigma(count).^(-2).*((kappa(count)-rho(count).*sigma(count)*1i.*u-d).*T(count)-2*log((1-g.*exp(-d.*T(count)))./(1-g)));
C = v0(count).*sigma(count).^(-2).*(kappa(count)-rho(count).*sigma(count)*1i.*u-d).*(1-exp(-d.*T(count)))./(1-g.*exp(-d.*T(count)));
charFunc = exp(A + B + C); 
FFTFunc = charFunc*exp(-rd(count)*T(count))./(alpha^2 + phi.*alpha - v.^2 + 1i*(phi.*2*alpha +1).*v);
payoff = real(exp(-1i.*v.*k).*FFTFunc);
end