The Merton Jump-Diffusion Model
Merton (1976)  jump-diffusion model is an extension to the Geometric Brownian Motion model, with the underlying asset exhibit jumps in addition to have continuous diffusion paths. The asset price evolves as
where μ is the drift rate, σ the volatility of S and dWt the increment of a Wiener process. The independent Poisson process dq has a value of 1 with probability λdt and 0 otherwise.
The following is a ThetaML implementation of the Merton Jump-Diffusion Model.
Merton jump diffusion Process simulation
Implementation as ThetaML
model JumpDiffusionProcess %Jump diffusion model for stock price movement %Example input parameters: %S0 = 100, r = 0.03, sigma = 0.2, lambda = 1 %nu = -0.1, delta = 0.05, t = 1, ht = 1/12 import S0 "Initial stock price" import r "Risk-free interest rate" import sigma "Stock price volatility" import lambda "Jump intensity" import nu "Mean of jumps" import delta "Std. of jumps" import t "Simulation time horizon" import ht "Discretization time step" export S "Simulated stock prices" S = S0 %initialize stock prices loop t/ht %loop t/ht times theta ht %theta advances time by ht units P = poissrnd(lambda*ht,size(rand())) U = exp(P*nu + sqrt(P)*delta*randn()) S = S * exp((r-lambda*(exp(nu+0.5*delta^2)-1)-0.5*sigma^2)*ht + sigma*sqrt(ht)*randn())*U end end
Note that this implementation requires the Matlab statistics toolbox. The model JumpDiffusionProcess uses Matlab function poissrnd to sample a poission distributed random variable.
Here is an implementation as a Theta.m function without using the Matlab statistics toolbox:
Implementation as Theta.m
function state = Theta(dt, state) if nargin == 0 state.S.value = 100; state.EUR.value = 1; else r = 0.05; mu1 = -0.2; lambda = 0.1; sigma = 0.2; gamma = 0.45; mu = r -(exp(mu1 + 0.5*gamma^2) - 1) * lambda; max_dt = lambda; n = max(dt/max_dt, 1); sub_dt = dt/n; dq = zeros(size(state.S)); for i = 1 : n jumpsize = mu1 + gamma*randn(size(state.S));% u = rand(size(state.S)); dq(u<lambda*sub_dt) = dq(u<lambda*sub_dt) + jumpsize(u<lambda*sub_dt); end state.S = state.S .* exp((mu-0.5*sigma^2)*dt + sqrt(dt)*sigma*randn(size(state.S))+dq); state.EUR = state.EUR* exp(-r*dt); end end
Analytic pricing of European calls
European call and put options with the underlying following a Merton jump-diffusion process have analytical solutions. Below is a Matlab implementation of the analytical solution for European calls.
function c = CPMertonJD % function c = CPMertonJD calculates the analytical price % of a European call under jump diffusion. So = 100; % Initial Stock price X = 100; % Strike T = 1; % Time to maturity delta = 0.4; % std of lognormal jump process nu = 0; % mean of lognormal jump process lambda = 5; % intensity of jumps vola = 0.4; % vola of stock price r = 0.05; % interest rate % The price is gven as a series of terms, % compute the first N terms for the result N = 50; K = exp(nu + 0.5*delta^2) - 1; c = 0; for n = 0 : N sigma_n = sqrt(vola^2 + n*delta^2/T); r_n = r - lambda*K + n*log(1+K)/T; d1 = (log(So/X) + (r_n + 0.5*sigma_n^2)*T)/(sigma_n*sqrt(T)); d2 = d1 - sigma_n*sqrt(T); f_n = So * normcdf(d1,0,1) - X * exp(-r_n*T) * normcdf(d2,0,1); c = c + exp(-lambda*(1+K)*T) * (lambda*(1+K)*T)^n * f_n/factorial(n); end end
Some observations about Monte Carlo convergence:
- The error decreases with increasing number of paths.
- The error is less if the variance of the jump-size is less.
- If the jump intensity (lambda) is high, then the error decreases at a smaller rate with increasing number of paths.
- ↑ Merton, R., Option pricing when underlying stock returns are discontinu- ous, Journal of Financial Economics, 3 (1976), p125–144.