Option pricing in Merton model, comparison between Merton series and Carr-Madan

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

%% 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
• Can you recover the analytical Black/Scholes solution when $\lambda = 0$? – LocalVolatility Nov 9 '17 at 17:11
• Not sure I understood your question, if $\lambda=0$ both the analytical solution (merton series) and the carr-madan solution are 0.132696765846609 (with other parameters as in the code above) – sound wave Nov 9 '17 at 23:23
• The background of my question is to help you narrow down where the error in your code comes from. The reason why I asked for $\lambda = 0$ is that you can't be 100% sure that your analytical pricer is correct and this reduces to the Black/Scholes case which you can independently check. Do some more checks - do the pricers still agree when $\lambda = 0$ but $K \neq S$ and so on. Eventually you should be able to narrow down the problem to a few candidate lines. – LocalVolatility Nov 9 '17 at 23:35
• Here are some reference values that should help you track down the errors. I am using the same parameters as you do in your code - all are calls for a fixed spot of S = 1. (1) K = 1, r = 10% -> C = 0.239266, (2) K = 1, r = 0% -> C = 0.202570, (3) K = 0.8, r = 10% -> C = 0.340487, (4) K = 0.8, r = 0% -> C = 0.292532, (5) K = 1.2, r = 10% -> C = 0.176873, (6) K = 1.2, r = 0% -> C = 0.150289. – LocalVolatility Nov 10 '17 at 0:00
• Thanks for the help, I got your same values (1) 0.239266363231017 (2) 0.202570314636424 (3) 0.340486886102500 (4) 0.292531648868869 (5) 0.176872747835715 (6) 0.150288774633102 – sound wave Nov 10 '17 at 3:32