In my latest book, I have a chapter titled “A Closer Look at the Square-root and 3/2-Model”.  One of the things I do in that chapter is give the joint transition density for the Heston ’93 model (which is associated to the Square-root model) and for the 3/2-model. These two models are popular stochastic volatility models (abbreviated ‘SV’).

A reader over at the Wilmott forums asked the natural question: what is the corresponding joint transition density when one adds price jumps to the Heston model? The resulting class of models is often abbreviated ‘SVJ’, standing for ‘stochastic volatility + jumps’.

As it turns out, the modifying the book formula for price jumps is quite easy. Showing that is the subject of this post and the next one. In this first one, I’ll just give the result. In the second one, I’ll give the derivation.

First, let’s review the results given in the book without the jumps.

The Joint Transition Density for the Heston Model

The Heston model is the SV system:

\displaystyle \left \{ \begin{array}{l}dS_t = b \, S_t \, dt + \sqrt{V_t} S_t \, dB_t, \\dV_t = (\omega - \theta V_t) \, dt + \xi \, \sqrt{V_t} \, dW_t, \quad dB_t \, dW_t = \rho \, dt. \end{array}\right. \quad\quad\quad (1)

Actually, we are considering two slightly different models at once: the “risk-neutral” version and the “real-world” version. In the risk-neutral model, b = r - q, where r is an interest rate and q is a dividend yield. Under the risk-neutral model, e^{-(r-q) t} \, S_t is a martingale. Under the real-world model, b is simply a constant parameter that one estimates. With those two possibilities for b, the joint transition density is given by

p(t,S_t,V_t|S_0,V_0) = \frac{1}{\pi S_t} \int_0^{\infty} \Re \left\{(\frac{S_t e^{-b t}}{S_0})^{-i u} G(t;-u,V_0,V_t) \right\} \, du, \quad (2)

where

G(t;z,v_0,v_t) = \frac{d \, e^{(1+\nu) \hat{\theta} \tau_t/2}}{2 \sinh(d \tau_t/2)} \, e^{-B_t v_0 - C_t v_t}  \left( \frac{v_t}{v_0} \right)^{\nu/2} I_{\nu} \left( \frac{d \sqrt{v_0 v_t}}{\sinh(d \tau_t/2)} \right),

using

\nu = \frac{2 \omega}{\xi^2} - 1, \quad \tau_t = \frac{1}{2}\xi^2 t, \quad B_t = \frac{r_1 - r_2 e^{d \tau_t}}{e^{d \tau_t}  - 1}, \quad C_t = \frac{r_1 e^{d \tau_t} - r_2}{e^{d \tau_t}  - 1}, 

\hat{\theta} = \frac{2}{\xi^2}(\theta + i z \rho \, \xi), \quad \tilde{c} = \frac{1}{\xi^2}(z^2 - i z), \quad d = (\hat{\theta}^2 + 4 \tilde{c})^{1/2}, \quad r_{1,2} = \frac{1}{2}(\hat{\theta} \pm d).

Equation (2) is equation (7.40) in the book. It looks messy, but essentially it’s just a single (numerical) integration involving elementary functions and the modified Bessel function: I_{\nu}(z). It’s quite straight-forward to implement in Mathematica, with the caveat that one should re-work the answer slightly to use only e^{-d \tau_t} terms, which helps prevent unexpected branch-cut crossings. The code, which uses that re-work, in given in the book in Sec. 7.12 Appendix.

In addition to ‘b’, note that the other parameters will generally be different under the two versions of the model: real-world or risk-neutral. The exception is \xi, which must be the same, as well as the entire process \{V_t\}, which also must be the same.

Applications

Once you have the joint density, you can compute various interesting things. For example, with the risk-neutral version you could compute the value of a call option on a stock with a single discrete dividend prior to expiration, now accounting for stochastic volatility effects. Or, you could compute the value of a compound option. If you’re still reading this, you probably can think of your own applications. With the real-world version, a time series of underlying security prices and a good proxy for the instantaneous volatility, you can make parameter estimates, using Maximum Likelihood Estimation.

The Joint Transition Density for the SVJ class of Models

For our purposes, we will define the ‘SVJ class of models’ to mean the Heston model modified to include independent, Poisson-driven, price jumps. Specifically, the modified SDE system becomes:

\displaystyle \left \{ \begin{array}{l}dS_t = b \, S_t \, dt + \sqrt{V_t} S_t \, dB_t + S_t (e^{\Delta X_t} - 1) dN_t, \\dV_t = (\omega - \theta V_t) \, dt + \xi \, \sqrt{V_t} \, dW_t, \quad dB_t \, dW_t = \rho \, dt. \end{array}\right. \quad\quad\quad (3)

The new notations include N_t, a Poisson process which counts the number of jumps that have occurred. Letting \lambda be the corresponding Poisson intensity, we have E_t[dN_t] = \lambda \, dt in SDE language, where E_t[\cdots] is the time t expectation, conditional on the state of the system (S_t,V_t). (Again, with all the parameters, there is a ‘real-world’ value and a `risk-neutral’ value, but we are suppressing that in a notational abuse)

In addition, given that there is a jump at time t, the jump size (actually the change in \log S_t) is denoted by \Delta X_t. This is a random variate with support on (-\infty,\infty). It is drawn independently from some jump density p_J(x) (real-world) or q_J(x) (risk-neutral). In a brief notation for this draw:  \Delta X_t \sim p_J(\cdot), and similarly for the risk-neutral case.

Our discussion allows p_J(\cdot)  or q_J(\cdot) to be rather general. What are the restrictions? First, the q_J(\cdot) density is subject to the condition that the integral in equation (4) below is well-defined.

Second , a  requirement from “no-arbitrage” considerations is that the two densities have the same support. For example, if there is a positive probability density that the stock market can crash (i.e., jump downward) by, let’s say, -20%, then there must also be a positive risk-neutral probability for the same event. However, the magnitudes can be, and usually are, quite different. Indeed, because investors are risk-averse, the risk-neutral probability of a broad-market crash is typically much higher than a plausible real-world probability estimate. In older literature, these risk-neutral probabilities are more aptly named ‘util-probs’, for ‘utility-adjusted probabilities’. See the mentioned book (Ch. 16, “Fear of Jumps”) for more on that utility connection.

Let’s note a special case. Namely, when both jump-size densities are normal densities, the resulting model is the Bates ’96 model, which blends the Heston model and Merton’s jump-diffusion model.

Another important point is that the constant ‘b’ in (3) has now changed under both interpretations: real-world or risk-neutral. That is, b is still a constant, but it’s numerical value is now generally going to be different from what it was under the pure Heston model. Under the real-world interpretation of the process, suppose you fix a parameterized density p_J(\cdot). Then, you go and estimate all the parameters by Maximum Likelihood Estimation. The estimate for b will change (even on the same observations) because the model has changed.

Under the risk-neutral interpretation of the process, we still need to enforce that S_t e^{-(r-q)t} is a martingale. Taking the time-t expectation of the first line of (3) yields the martingale condition:

\displaystyle r - q = b + \lambda \int_{-\infty}^{\infty} (e^x - 1) \, q_J(x) \, dx = b + \lambda (\hat{q}_J(-i) -1). \quad\quad (4)

Note that I have introduced the Fourier transform of the (risk-neutral) jump-size density: \hat{q}_J(z) = \int_{-\infty}^{\infty} e^{i z x} \, q_J(x) \, dx. In a similar notation, \hat{p}_J(z) denotes the transform of the real-world density.

Finally, using these notations, what is the modification of (2) that incorporates jumps? It turns out to be quite simple: just insert a relatively simple factor into the integrand:

p(t,S_t,V_t|S_0,V_0) = \frac{1}{\pi S_t} \int_0^{\infty} \Re \left\{(\frac{S_t e^{-b t}}{S_0})^{-i u} \, e^{\psi(u) t} \, G(t;-u,V_0,V_t) \right\} \, du, \quad (5)

where

\displaystyle \psi(u) = \left \{ \begin{array}{ll} \lambda  \, (\hat{p}_J(u) - 1), & \mbox{(real-world)}, \\  -i u \lambda \, (\hat{q}_J(-i) - 1) + \lambda \, (\hat{q}_J(u) - 1), & \mbox{(risk-neutral)}. \end{array}\right. 

In the risk-neutral case, you keep b = r - q in (5). Then the additional term in \psi(u) accounts for the modification of b given in equation (4).

For the special case of the Bates model, we can write the jump term in (3) as

e^{\Delta X_t} - 1 = e^{\alpha + \delta \epsilon_t} -1,

where \epsilon_t is a draw from an independent standard normal distribution. Thus the (log) jump-sizes are normally distributed with mean \alpha and variance \delta^2. In that case (5) uses

\displaystyle \psi(u) = \left \{ \begin{array}{ll} \lambda  \, (e^{i u \alpha - u^2 \delta^2/2} - 1), & \mbox{(real-world)}, \\ -i u \lambda \, (e^{ \alpha + \delta^2/2} - 1) + \lambda \, (e^{i u \alpha - u^2 \delta^2/2} - 1), & \mbox{(risk-neutral)}. \end{array}\right. 

Note that the values of \lambda, \alpha, \, \mbox{and} \, \delta will general differ between the real-world and risk-neutral cases, although my notation does not explicitly indicate this.

As it turns out, the modification here to add price jumps to the risk-neutral model is exactly the same as the modification discussed in J. Gatheral’s book “The Volatility Surface” (pg. 66) for the characteristic function that prices vanilla options. This is not too surprising — but deserves to be shown carefully: perhaps I will do that in a future post.