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);
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
  • $\begingroup$ Can you recover the analytical Black/Scholes solution when $\lambda = 0$? $\endgroup$ Commented Nov 9, 2017 at 17:11
  • $\begingroup$ 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) $\endgroup$
    – sound wave
    Commented Nov 9, 2017 at 23:23
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Nov 9, 2017 at 23:35
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Nov 10, 2017 at 0:00
  • $\begingroup$ 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 $\endgroup$
    – sound wave
    Commented Nov 10, 2017 at 3:32


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.