The Merton Jump Diffusion Model

Merton Jump Diffusion Model

The Merton Jump Diffusion model proposed by Merton in 1976 is an extension of the Black-Scholes model (link).

It contains:

• A diffusion part with a standard wiener process W(t) (link):

μ: annualised expected return on the asset price, σ: annualised volatility of the asset price, W_t: Brownian motion.

• A jump part with a compound Poisson process Q(t):

Q_t: compound Poisson process

Y_i: price ratio associated with the i-th jump happening at the random time t_i. (Y_i) are i.i.d, following a lognormal distribution, and are independent of W_t and N_t, N_t being a Poisson process (link) with intensity λ. E(Q_t) = λ.k.t, so Q_t – λ.k.t is a martingale.

Reference article: Merton, R. C. (1976) “Option Pricing When underlying stock returns are discontinuous”, Journal of Financial Economics, 3, 125-144

Compound Poisson Process

A compound Poisson process is a stochastic process with jumps. The arrival of events (jumps) follows a Poisson process (link). In addition to the random arrival time of events, the magnitude of the jumps is also random, with a specified probability distribution. The magnitudes of the events are independent of each other and also independent of the arrival times.

Example:

λ = 10, T = 1y, jumps ~ N(0, 0.2²). 10 jumps per year on average, the magnitude of the jumps follows a normal distribution N(0, 0.2²)

Merton Jump Diffusion Model: Explicit Solution

In the Merton Jump Diffusion model, the asset price S_t follows the SDE:

There is an explicit solution of this SDE¹:

¹Merton, R. C. (1976) “Option Pricing When underlying stock returns are discontinuous”, Journal of Financial Economics, 3, 125-144

Merton Jump Diffusion Model: Simulation

We can simulate the asset price:

S₀ = 100, μ = 5%, σ = 20%, λ = 10, T = 5y, jumps ~ N(0, 0.1²), δt = 5 / 1000

Probability Density Function of the Log Returns

If there are n jumps between 0 and t, we have: 

Knowing the number of jumps we can easily deduce the density of ln(S_t/S₀) which follows a normal distribution:

And we obtain the density of ln(S_t/S₀) with the following series using the probability to have n jumps and the conditional distribution:

How Jump Parameters Impact the Moments of the Distribution

The probability density function of the log-returns is not normal when there are jumps. In addition to an increase of the standard deviation of the log-returns, adding jumps adds tail-risk (kurtosis) and asymmetry (skewness), positive or negative depending on the sign of the average jump.

Jumps	Distribution jumps	St Deviation	Skewness	Excess Kurtosis
No (λ = 0)		+	0	0
Yes (λ > 0)	m = 0, δ > 0	++	0	+
Yes (λ > 0)	m > 0, δ > 0	++	+	+
Yes (λ > 0)	m < 0, δ > 0	++	–	+

μ = 5%, σ = 20%, λ = 0 or 1, T = 1y, m = 0 or -0.2, δ = 0.2

Save 10% on All Quant Next Courses with the Coupon Code: QuantNextBlog10

For students and graduates: We offer a 50% discount on all courses, please contact us if you are interested: contact@quant-next.com

Python Code

Import Libraries

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
from scipy.stats import norm
import math

Compound Poisson Process

#Simulate Normal Compound Poisson Process
def normal_compound_poisson_process(lambda_, T, m, delta):

    t = 0
    jumps = 0
    event_values = []
    event_times = []
    event_values.append(0)
    event_times.append(0)
    while t < T:
        t = t + np.random.exponential(1 / lambda_)
        jumps = jumps + np.random.normal(m, delta)
        if t < T:
            event_values.append(jumps)
            event_times.append(t)

    return event_values, event_times

#Plot Normal Compound Poisson Process
def plot_normal_compound_poisson_process(lambda_, T, m, delta):
    event_values, event_times = normal_compound_poisson_process(lambda_, T, m, delta)

    plt.figure(figsize=(10, 8))
    plt.step(event_times, event_values, where='post')
    plt.xlabel('Time (in years)')
    plt.ylabel('Value')
    plt.title('Normal Compound Poisson Process Simulation')
    plt.grid(True)
    plt.show()

#Example
lambda_ = 10  # Intensity of Poisson process
T = 1          # Time
m = .0       # Mean of jumps
delta = .2    # Standard deviation of jumps

plot_normal_compound_poisson_process(lambda_, T, m, delta)

Merton Jump Diffusion Process

#Simulate Merton Jump Diffusion (MJD) Process
def MJD_process(S0, mu, sigma, lambda_, T, m, delta, num_steps):

    dT = T / num_steps
    price = np.zeros(num_steps + 1)
    price[0] = S0
    k = np.exp(m + .5 *delta**2) - 1
    
    for t in range(1, num_steps + 1):
        # Brownian motion
        Z = np.random.normal(0, 1)
        dW = np.sqrt(dt) * Z
        #jumps between t and t + dt
        jumps = np.sum(normal_compound_poisson_process(lambda_, dT, m, delta)[0])
        price[t] = price[t - 1] * np.exp((r - .5 * sigma**2 - lambda_ * k) * dt + sigma * dW + jumps)
        
    return price

#Plot Merton Jump Diffusion Process
def plot_MJD_process(S0, mu, sigma, lambda_, T, m, delta, num_steps):

    price_simu = MJD_process(S0, mu, sigma, lambda_, T, m, delta, num_steps)
    time_axis = np.linspace(0, T, num_steps + 1)

    plt.figure(figsize=(10, 5))
    plt.step(time_axis, price_simu, where='post')
    plt.xlabel('Time (in years)')
    plt.ylabel('Price')
    plt.title('Merton Jump Diffusion Process Simulation')
    plt.grid(True)
    plt.show()

#Example
S0 = 100   #initial price
mu = .05   #expected return
sigma = .2   #volatility
lambda_ = 10  # Intensity of Poisson process
T = 5          # Time
m = 0.0       # Mean of jumps
delta = .1    # Standard deviation of jumps
num_steps = 1000 #number of time steps

plot_MJD_process(S0, mu, sigma, lambda_, T, m, delta, num_steps)

Probability Density Function Log Return with Merton Jump Diffusion Model

#Density Log Return Merton Jump Diffusion (MJD)
def MJD_density(x, mu, sigma, lambda_, T, m, delta, nmax):
    res = 0
    k = np.exp(m + .5 *delta**2) - 1
    for n in range(nmax):
        proba_n = np.exp(-lambda_ * T) * (lambda_ * T)**n / math.factorial(n) #proba n jumps
        a_n = (mu - .5 * sigma**2 - lambda_ * k) * T + n * m #avg conditional normal distribution
        s_n = (sigma**2 * T + n * delta**2)**.5
        res = res + proba_n * norm.pdf(x, a_n, s_n)
    return res
                                       
#Plot Density Log Return MJD
def plot_MJD_density(mu, sigma, lambda_, T, m, delta, nmax):
    x = np.linspace(-3., 3., 100)
    params = [mu, sigma, lambda_, T, m, delta, nmax]
    density_x = [MJD_density(x_, *params) for x_ in x]

    plt.figure(figsize=(10, 5))
    plt.plot(x, density_x)
    plt.xlabel('Log Return')
    plt.ylabel('Density')
    plt.title('Density Log Return with Merton Jump Diffusion')
    plt.grid(True)
    plt.show()
    

#Example
mu = .05   #expected return
sigma = .2   #volatility
lambda_ = 10  # Intensity of Poisson process
T = 1          # Time
m = 0.0       # Mean of jumps
delta = .1    # Standard deviation of jumps
nmax = 100 #maximum number of jumps

plot_MJD_density(mu, sigma, lambda_, T, m, delta, nmax)

#Example
mu = .05   #expected return
sigma = .2   #volatility
#lambda_ = 10  # Intensity of Poisson process
T = 1          # Time
m = 0.0       # Mean of jumps
delta = .2    # Standard deviation of jumps
nmax = 100 #maximum number of jumps

x = np.linspace(-3., 3., 100)
params0 = [mu, sigma, 0, T, m, delta, nmax]
params1 = [mu, sigma, 1, T, m, delta, nmax]
params2 = [mu, sigma, 1, T, -0.2, delta, nmax]
density_x0 = [MJD_density(x_, *params0) for x_ in x]
density_x1 = [MJD_density(x_, *params1) for x_ in x]
density_x2 = [MJD_density(x_, *params2) for x_ in x]
plt.figure(figsize=(10, 5))
plt.plot(x, density_x0, label = 'Normal Distribution (No Jumps)')
plt.plot(x, density_x1, label = 'Distribution with Jumps')
plt.plot(x, density_x2, label = 'Distribution with Jumps, m < 0')
plt.xlabel('Log Return')
plt.ylabel('Density')
plt.title('Density Log Return with Merton Jump Diffusion')
plt.grid(True)
plt.legend()
plt.show()