I'm studying the Merton model for pricing an European call option. The jump-diffusion process is: $$X_t=bt+\sigma W_t+\displaystyle\sum_{i=1}^{N_t}Y_i.$$
$N_t$ is the Poisson process,
$W_t$ is the brownian motion with drift $b=\left(r-\frac{\sigma^2}{2}-\lambda\left(e^{\mu+\delta^2/2}-1\right)\right)$,
$\mu$ is the mean of jumps,
$\delta$ is the standard deviation of jumps,
$\Phi_X(u) = exp\left(t\left(ibu-\frac{\sigma^2u^2}{2}+\lambda\left(exp\left(i\mu u-\frac{\delta^2u^2}{2}\right)-1\right)\right)\right)$ is the characteristic function of $X_t$,
$V_0 = \displaystyle\sum_{n=0}^{+\infty}e^{-\lambda T}\frac{(\lambda T)^n}{n!}V_{BS}$ is the Merton series for the option price, where $V_{BS}=blsprice$ is the option price in Black-Scholes. So $V_0$ is the analytical solution.
I already solved the PIDE, derived from the jump-diffusion model, with an implicit scheme. The numerical solution approximates very well the analytical one.
Now I'm studying the Carr-Madan algorithm on the article by Tankov and Voltchkovato (pages 9 to 11), and I would like to compare its numerical solution with the analytical one of the Merton series.
Carr-Madan is based on the idea of considering the log of the strike, $k=log(K)$, and then compute the Fourier transform of the difference between the price of an European call in the jump-diffusion model (Merton model here)(this price is what we want to achieve) and the price of an European call in the Black-Scholes model. This Fourier transform is given by $$\zeta_T(v)=e^{ivrT}\frac{\Phi_X(v-i)-\Phi_{BS}(v-i)}{iv(1+iv)}$$ where $\Phi_X$ is the one above and $\Phi_{BS}=exp\left(-\frac{\Sigma^2T}{2}\left(v^2+iv\right)\right)$ is the characteristic function of log-stock price in the Black-Scholes model (in my code I'm using $\Sigma=\sigma$).
The call price in the Merton model can be computed by evaluating numerically the inverse Fourier transform of $\zeta_T$ and then by summing to it the Black-Scholes call price.
The inverse Fourier transform of $\zeta_T$ is (truncate and discretize the integral) $$\frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-ikv}\zeta_T(v)dv \approx \frac{L}{2\pi(N-1)} \sum_{m=0}^{N-1} w_m\zeta_T(v_m) e^{-ikv_m}$$ where $v_m = -L/2+m\Delta,~\Delta = L/(N-1)$ is the discretization step and $w_m$ are weights for a integration rule (I'm using Simpson's rule).
We can identify this sum with the DFT $F_n = \sum_{k=0}^{N-1} f_k e^{-2\pi ink/N}$ by chosing the log strikes $k_n = k_0+2\pi n/(N\Delta),~n=0,...,N-1$: $$\frac{L}{2\pi(N-1)} e^{ik_nL/2} \sum_{m=0}^{N-1} \underbrace{w_m\zeta_T(k_m) e^{-ik_0m\Delta}}_\textrm{$f_m$} e^{-2\pi inm/N}.$$
I'm studying the particular case with $\mu=0$ and $K=1$, that's why I compare (in the last line of the code) the analyrical solution with the first component of the Carr-Madan solution (first strike for Carr-Madan is $K=e^{k_0}=1$ because $k_0=0$).
But I'm facing some problems:
(1) if the initial underlying price $S_0=1$, then Carr-Madan solution is about $8$% bigger than the analytical one;
(2) if $S_0=1$ AND I remove the risk-free interest $r$ from $b$, then the Carr-Madan solution is equal to the analytical one (first seven digits after the separator are equal);
(3) if $S_0\neq 1$ the Carr-Madan solution is different from the analytical one (the error is smaller if I eliminate $r$ from $b$).
I really don't understand this behaviour and where I made the mistakes in the code, so I post it here. For the variables I used the names as in Tankov-Voltchkovato article. So this is the matlab code:
clear all
close all
format long
%% parameters
t = 0; % now
T = 1; % expiry date
tau = T - t; % time to expiry
S_0 = 1; % initial underlying price
r = 0.1; % risk-free rate
sigma = 0.2; % volatility
lambda = 1; % mean no. of jumps per unit time
mu = 0; % mean of jumps
delta = 0.5; % standard deviation of jumps
%% merton series (analytical solution)
K = 1; % strike
N = 30;
v = zeros(N+1,1);
for n = 0:N
sigma_n = sqrt(sigma^2+n*delta^2/tau);
Sn = S_0*exp(n*mu+n*delta^2/2-lambda*tau*exp(mu+delta^2/2)+lambda*tau);
v(n+1) = exp(-lambda*tau)*(lambda*tau)^n/factorial(n)*blsprice(Sn,K,r,tau,sigma_n);
end
MertonSeriesValue = sum(v);
%% carr-madan (numerical solution)
L = 2^16; % truncation interval
N = 2^16; % no. of points
v_m = linspace(-L/2,L/2,N);
Delta = L/(N-1); % discretization step
% corrected Simpson-Cavalieri weight rule
w = 4*ones(1,N);
w(1:2:N) = 2;
w(1) = 1;
w = w/3;
% log strikes
k_0 = 0;
k_n = k_0+(2*pi*(0:(N-1)))/(Delta*N);
K_CR = exp(k_n);
zeta_T = fourier(v_m,T,sigma,r,lambda,mu,delta); % Fourier transform
dft = fft(w.*zeta_T.*exp(-1i*k_0*Delta*(0:(N-1)))); % FFT computes DFT
z_T = L/(2*pi*(N-1))*exp(1i*k_n*L/2).*dft; % aproximate option prices
C_BS = blsprice(S_0,K_CR,r,T,sigma); % option value for Black-Scholes
C = z_T+C_BS; % option value for Merton
CarrMadanMertonPrice = max(real(C),0);
%% comparison
Merton_CR_KK = [MertonSeriesValue(1); CarrMadanMertonPrice(1); K_CR(1)==K]
function used at line 46
function[zeta_T] = fourier(v,T,sigma,r,lambda,mu,delta)
%
% Fourier transform in log-strike k of z_T(k)
%
u = v-1i;
b = r-(sigma^2)/2+lambda*(1-exp(mu+(delta^2)/2)); % if I remove r from here the solution will be more precise
phi_T_sigma = exp(-(sigma^2)*T*(u.^2+1i*u)/2); % characteristic function of log-stock price in the Black-Scholes model
phi_T = exp(T*(-((sigma*u).^2)/2+1i*b*u+lambda*(exp(-(delta*u).^2/2+1i*mu*u)-1))); % characteristic function of the jump-diffusion process in the merton model
zeta_T = exp(1i*v*r*T).*(phi_T-phi_T_sigma)./(1i*v.*(1+1i*v)); % Fourier transform
