All posts by Alan

Errata for ‘Back to Basics: a new approach to the discrete dividend problem’

Here is a small correction to Table 6 in this 2003 article by E.G. Haug, J. Haug, and A. Lewis.

(Published value)
(Corrected value)

Table 6 shows various approximations to the value of a European call option under discrete dividends, with multiple dividends. The table column labeled ‘HHL’ is referred to as an ‘exact integration’. This language suggests that all the displayed digits should be correct. However, a careful reader was attempting to replicate those values and did not agree with all the digits. (Thank you, Yiannis!). Rerunning the codes a little more carefully yields new best estimates shown above. Note the T = 1 entry is unchanged. The revised entries also agree, to the digits shown, with the mentioned reader’s independent calculations, which use a different method.

I frequently refer to this article in Ch. 9 of ‘Option Valuation under Stochastic Volatility II’. In particular, under ‘Suggested reader projects’ I suggest the reader extend the chapter numerics to call options and compare results to the HHL tables. So, for those who might be doing that, the above small revisions may be meaningful.


Espen G. Haug, Jørgen Haug, and Alan Lewis. Back to Basics: a new approach to the discrete dividend problem. Wilmott magazine, pgs. 37-47, Sept. 2003

The Joint Transition Density for SVJ Models — Part I

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’).

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)


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),


 \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.


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)


\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: see the next post for that.

A Young Quant’s Illustrated Primer: the SABR Model

What is the SABR model?

The SABR model (pronounced as in light saber) is a popular stochastic volatility model. A somewhat more general version than the one we discuss here was introduced by P. Hagan and co-authors in 2002, with applications to interest rate-related derivatives. The model became popular because a (quite difficult) small-time analysis resulted in some easy-to-apply, approximate, “smile” formulas that could fit the authors’ target markets.

While equity applications are rare, the model makes mathematical sense for equities, where the stochastic process pair is \{S_t,\sigma_t\}. Here  S_t \ge 0 is a stock price and \sigma_t > 0 is the associated stochastic volatility. Using that equity language, the SABR model, as an SDE (stochastic differential equation) system, reads

\displaystyle \left \{ \begin{array}{l}dS_t = \sigma_t S_t^{\beta} dB_t, \quad (t < \tau_0), \\d \sigma_t = \nu \, \sigma_t \, dW_t, \quad dB_t \, dW_t = \rho \, dt. \end{array}\right. \quad\quad\quad\quad\quad (1)

In equation (1), \beta is any real number, and \nu > 0 is the “volatility of volatility”. The latter can be set equal to one without loss of generality, so we’ll do so for the remainder of this post. Although any real \beta makes mathematical sense, fits to market data will almost always have -\infty < \beta \le 1, so we’ll stick to that range also. Also \rho \in (-1,1) is a correlation parameter. Finally, \tau_0 is the first time, if any, that the stock price hits 0.

Boundary conditions

The stock price can hit 0 with a positive probability whenever \beta is strictly less than 1. What happens after that requires a boundary condition. Of course, for equities, to avoid an arbitrage opportunity, once a stock price hits zero, it must stay there (unless the company is recapitalized). So, an additional specification of the equity version of the model is that S_t = 0 for all t \ge \tau_0.

For the original fixed income applications, the underlying security is a (forward) interest rate and other boundary conditions are possible.

Exact solutions

By a “solution” to equation (1), I mean a solution for the associated transition probability density. Once you have that, option values, for example, are found by an integration. Exact solutions are known for the following special cases:

  • Normal SABR model: \beta = 0 with \rho = 0.
  • Lognormal SABR model: \beta = 1 with any \rho.

McKean’s formula

The Normal SABR model solution is found by adapting a famous formula for the transition density found by Henry P. McKean in 1970. For that case, it makes more sense to switch gears slightly in the notation and use X_t instead of S_t, now X_t the location of a “particle” which can hit 0 and keep on going to arbitrarily negative values. With that new notation, the Normal SABR model process is:

\displaystyle \left \{ \begin{array}{l}dX_t = \sigma_t dB_t, \\d \sigma_t = \sigma_t \, dW_t. \end{array}\right. \quad\quad\quad\quad (2)

Now B_t and W_t are two independent Brownian motions. Equation (2) also describes a generalized type of Brownian motion, namely Brownian motion on a two-dimensional hyperbolic space.

Suppose we start the particle at X_0 = x_0 = 0 and \sigma_0 = 1. Then, the probability of finding the particle at subsequent times is given by Mckean’s formula for G(t,X,\sigma|x_0,\sigma_0). If we let t range from 0.01 to 1, we can plot McKean’s formula to find the following evolution:

Evolution of McKean formula probability density for the Normal SABR model
Evolution of McKean formula probability density for the Normal SABR model

If you watch closely, you’ll see that the probability density starts with a spike near (X,\sigma) = (0,1) and then spreads out in the X direction. In the \sigma direction, mass begins to accumulate near \sigma \rightarrow 0. That’s due to lack of mean reversion in the volatility process and is an ‘unphysical’ aspect of the model.

Now that we’ve talked about it so much, you may be curious as to what this mysterious formula looks like. Here it is:

\displaystyle G(t,X,Y|x,y) = \frac{e^{-t/8}}{2 \, Y^2 \, (\pi t)^{3/2}} \int_d^{\infty} \frac{s \, e^{-s^2/(2 t)}}{\sqrt{\cosh s - \cosh d}} \, ds,

\displaystyle \mbox{where} \quad d = \cosh^{-1}z, \quad \mbox{using} \quad z = 1 + \frac{(X-x)^2 + (Y-y)^2}{2 y Y}

To adapt the McKean solution for an equity model, you can use the well-known “method of images”. The method will ensure that the associated equity process S_t, which starts at S_0 > 0, is absorbed should it hit S=0.

(There’s much more to the model than we have touched upon, and we’ll certainly return to it in future blog posts. This post has been adapted from material in Ch. 8 (“A Closer Look at the SABR Model”) from “Option Valuation under Stochastic Volatility II”).

Further reading and codes

Animation code

Started adding book codes from ‘Option Valuation under Stochastic Volatility II’

I have started adding book codes to this site. Almost all will be Mathematica notebooks (ending in .nb) , although there will be some C/C++. This may be a slow process since most files will require some clean-up before posting.

One file has been added so far: from Ch. 1: CIRJumpStart.nb, which was used to generate Fig. 1.6 in the book and below:

Fig. 1.6 CIR Bond Process with a Delayed Jump-start.
Fig. 1.6. CIR Bond Process with a Delayed Jump-start.

The figure shows a fit of the CIR bond process with a delayed jump-start (solid curve) to the US Treasury Zero Coupon yield curve (dashed curve), using Feb. 10, 2012 data. The delayed jump-start is/was a mechanism to model the end of the ZIRP (Zero Interest Rate Policy) regime. Details and related models are found in the book Chapter 1: “Slow Reflection, Jump-returns, and Short-term Interest Rates”.

The code file location is the Latest Book link (Just click on ‘Selected Codes’)

I will probably just work my way slowly through the book, but if you are a reader and particularly want to see certain book codes, please let me know by a comment or email. I will try to accommodate if possible.