As with the previous article on the Heston Stochastic Volatility model, this article will deal with another assumption of the Black-Scholes model and how to improve it. In the Black-Scholes model the stock price evolves as a geometric Brownian motion. Crucially, this allows continuous Delta hedging and thus a fixed no-arbitrage price for any option on the stock.

If we relax the assumption of GBM and introduce the concept of discontinuous jumps in the stock price then it becomes impossible to perfectly hedge and leads to *market incompleteness*. This means that options prices are only *bounded* rather than fixed. In this article we will consider the effect on options prices when such jumps occur and implement a semi-closed form pricer in C++ based on the analytical formula derived by Merton^{[2]}.

## Modelling Jump-Diffusion Processes

*This section closely follows the chapter on Jump Diffusions in Joshi ^{[1]}, where more theoretical details are provided.*

In order to model such stock "jumps" we require certain properties. Firstly, the jumps should occur in an instantaneous fashion, neglecting the possibility of a Delta hedge. Secondly, we require that the probability of any jump occuring in a particular interval of time should be approximately proportional to the length of that time interval. The statistical method that models such a situation is given by the Poisson process^{[1]}.

A Poisson process states the probability of an event occuring in a given time interval $\Delta t$ is given by $\lambda \Delta t + \epsilon$, where $\lambda$ is the *intensity* of the process and $\epsilon$ is an error term. The integer-valued number of events that have occured at time $t$ is given by $N(t)$. A necessary property is that the probability of a jump occuring is independent of the number of jumps already occured, i.e. all future jumps should have no "memory" of past jumps. The probability of $j$ jumps occuring by time $t$ is given by:

Thus, we simply need to modify our underlying GBM model for the stock via the addition of jumps. This is achieved by allowing the stock price to be multiplied by a random factor $J$:

\begin{eqnarray} d S_t = \mu S_t dt + \sigma S_t dW_t + (J-1) S_t dN(t) \end{eqnarray}Where dN(t) is Poisson distributed with factor $\lambda dt$.

## Prices of European Options under Jump-Diffusions

After an appropriate application of risk neutrality^{[1]} we have that $\log S_T$, the log of the final price of the stock at option expiry, is given by:

This is all we need to price European options under a Monte Carlo framework. To carry this out we simply generate multiple final spot prices by drawing from a normal distribution and a Poisson distribution, and then selecting the $J_j$ values to form the jumps. However, this is somewhat unsatisfactory as we are specifically choosing the $J_j$ values for the jumps. Shouldn't they themselves also be random variables distributed in some manner?

In 1976, Robert Merton^{[2]} was able to derive a semi-closed form solution for the price of European options where the jump values are themselves normally distributed. If the price of an option priced under Black-Scholes is given by $BS(S_0,\sigma,r,T,K)$ with $S_0$ initial spot, $\sigma$ constant volatility, $r$ constant risk-free rate, $T$ time to maturity and $K$ strike price, then in the jump-diffusion framework the price is given by^{[1]}:

Where

\begin{eqnarray} \sigma_n &=& \sqrt{\sigma^2 + n \nu^2 / T} \\ r_n &=& r - \lambda(m - 1) + n \log m / T \\ \lambda ' &=& \lambda m \end{eqnarray}The extra parameters $\nu$ and $m$ represent the standard deviation of the lognormal jump process and the scale factor for jump intensity, respectively.

## C++ Implementation

We will avoid a fully-fledged object-oriented approach for this model as it extends simply the analytical pricing of European call options quite straightforwardly. Below is the code in its entirety. The major change includes the calculation of the factorial within the summation loop and the weighted sum of the Black-Scholes prices:

```
#define _USE_MATH_DEFINES
#include <iostream>
#include <cmath>
// Standard normal probability density function
double norm_pdf(const double x) {
return (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x);
}
// An approximation to the cumulative distribution function
// for the standard normal distribution
// Note: This is a recursive function
double norm_cdf(const double x) {
double k = 1.0/(1.0 + 0.2316419*x);
double k_sum = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937 + k*(-1.821255978 + 1.330274429*k))));
if (x >= 0.0) {
return (1.0 - (1.0/(pow(2*M_PI,0.5)))*exp(-0.5*x*x) * k_sum);
} else {
return 1.0 - norm_cdf(-x);
}
}
// This calculates d_j, for j in {1,2}. This term appears in the closed
// form solution for the European call or put price
double d_j(const int j, const double S, const double K, const double r, const double v, const double T) {
return (log(S/K) + (r + (pow(-1,j-1))*0.5*v*v)*T)/(v*(pow(T,0.5)));
}
// Calculate the European vanilla call price based on
// underlying S, strike K, risk-free rate r, volatility of
// underlying sigma and time to maturity T
double bs_call_price(const double S, const double K, const double r,
const double sigma, const double T) {
return S * norm_cdf(d_j(1, S, K, r, sigma, T))-K*exp(-r*T) *
norm_cdf(d_j(2, S, K, r, sigma, T));
}
// Calculate the Merton jump-diffusion price based on
// a finite sum approximation to the infinite series
// solution, making use of the BS call price.
double bs_jd_call_price(const double S, const double K, const double r,
const double sigma, const double T, const int N, const double m,
const double lambda, const double nu) {
double price = 0.0; // Stores the final call price
double factorial = 1.0;
// Pre-calculate as much as possible
double lambda_p = lambda * m;
double lambda_p_T = lambda_p * T;
// Calculate the finite sum over N terms
for (int n=0; n<N; n++) {
double sigma_n = sqrt(sigma*sigma + n*nu*nu/T);
double r_n = r - lambda*(m - 1) + n*log(m)/T;
// Calculate n!
if (n == 0) {
factorial *= 1;
} else {
factorial *= n;
}
// Refine the jump price over the loop
price += ((exp(-lambda_p_T) * pow(lambda_p_T,n))/factorial) *
bs_call_price(S, K, r_n, sigma_n, T);
}
return price;
}
int main(int argc, char **argv) {
// First we create the parameter list
double S = 100.0; // Option price
double K = 100.0; // Strike price
double r = 0.05; // Risk-free rate (5%)
double v = 0.2; // Volatility of the underlying (20%)
double T = 1.0; // One year until expiry
int N = 50; // Terms in the finite sum approximation
double m = 1.083287; // Scale factor for J
double lambda = 1.0; // Intensity of jumps
double nu = 0.4; // Stdev of lognormal jump process
// Then we calculate the call jump-diffusion value
double call_jd = bs_jd_call_price(S, K, r, v, T, N, m, lambda, nu);
std::cout << "Call Price under JD: " << call_jd << std::endl;
return 0;
}
```

The output from the code is:

`Call Price under JD: 18.7336`

We can clearly see that in comparison to the Black-Scholes price of $10.4506$ given in this article, the value of the call under the jump diffusion process is much higher. This is to be expected since the jumps introduce extra volatility into the model.

In the next article we will consider the pricing of exotic options in a fully object-oriented environment using jump diffusion models.