Local stochastic volatility: shaken, not stirred

Dominique Bang introduces a novel LSV approach to term distribution modelling

CLICK HERE TO VIEW THE PDF

CLICK HERE TO LISTEN TO DOMINIQUE IN CONVERSATION WITH RISK.NET

For the past two decades, local stochastic volatility (LSV) models have been enjoying a great deal of popularity among practitioners, not only as term structure models for the pricing and risk management of complex products (see, for example, Andersen & Piterbarg 2010; Blacher 2001; Lipton 2002), but also for the modelling of vanilla options (see, for example, Berestycki et al (2004) and the widely used stochastic alpha, beta, rho (SABR) formula from Hagan et al (2002)). Non-LSV models have also attracted some interest, eg, the Lévy process class (see Ticot & Charvet (2013) for an application to inflation modelling), but are used far less in practical applications.

In the world of vanilla interest rate derivatives (a particular focus of this article), standard LSV frameworks have recently been put under increasing pressure from changing market environments, be it low or negative rates in major currencies or regulatory stress tests scenarios. To compound these difficulties, there are ever-growing demands on these models’ ability to calibrate to a large number of swaptions, while simultaneously capturing the convexity embedded in constant maturity swap (CMS) products accurately. As a result, the near-standard SABR model volatility marking approach, based on the small-time asymptotics in Hagan et al (2002), is now problematic, as the model is both insufficiently rich to match the market and also plagued by well-known arbitrage problems for high and low option strikes. A number of authors have proposed different approaches to tackle the current challenges. For example, Andreasen & Huge (2013) use an extended SABR model based on a one-step partial differential equation (PDE) method, producing an arbitrage-free interpolator, but with a term distribution that deviates from the original dynamics as time to maturity increases. Balland & Tran (2013) propose a semi-analytical approach that is not guaranteed to be arbitrage free (though it is stated that it will be in practical applications). Antonov et al (2017) present a mixture approach.

While the improvements to the original SABR framework represent meaningful progress, none is, in our opinion, completely satisfactory from a theoretical or a practical perspective, and many involve non-intuitive parametrisation, a possibility of arbitrage, a lack of flexibility or challenges in the implementation. Here, we present a new approach that allows us to mix a generic local volatility (LV) function, parametric or non-parametric, with a pure stochastic volatility (SV) process in an intuitive, tractable, easy-to-implement and arbitrage-free way, in the sense that there are no negative butterflies for a given maturity. One design strength of the approach is that the choice of the LV component can be guided by criteria based on relevance rather than numerical and analytical feasibility.

The rest of this article is organised as follows. We first describe the LSV dynamics, present the fundamental idea of our method and show how to handle the LV component (via a one-dimensional integral), provided that the pure SV process is tractable, either via its cumulative density function (CDF) or, equivalently, via its European options prices. Then, we illustrate the method, pairing a general LV with a normal SABR process. The normal SABR dynamics are somewhat tractable, allowing for option prices through a one-dimensional integral of special functions (see Antonov et al 2017); unfortunately, the existing formulas become computationally too expensive when embedded in our method and combined with an LV process. Thus, we develop new analytical formulas for option prices and probability densities in the normal SABR model. These approximations are arbitrage free by construction, computationally very simple (only a few calls to the Gaussian CDF are required) and accurate, even for maturities as long as 50 years. Next, we make a concrete choice for the LV function and eventually show that the model (which supports negative rates) can be successfully calibrated to a large strip of options and to CMS forward levels. Last, we consider, for variety, an (arbitrary) piecewise linear LV combined with a normal Heston, describe the calibration strategy and provide option-pricing formulas, as a discrete sum of options under the Heston model. Finally, we conclude the paper.

Local stochastic volatility

Notation and dynamics

Our primary objective in this article is to compute European option pricing on some forward process $F_t$, in the context of a general LSV model of the form:

	$\mathrm{d}{F}_t$	$=\sigma (F_t)\mathrm{d}{M}_t$		(1)
	$\mathrm{d}{M}_t$	$={\alpha }_t\mathrm{d}{W}_t$		(2)

In (1) and (2), the LV function $\sigma$, assumed positive, controls the support of the distribution and the dynamics of the model (also known as ‘backbone’). The pure SV component, $M_t$, is driven by a Brownian motion $W_t$, with some generic volatility process ${\alpha }_t$ that will typically depend on a second source of randomness. In most model configurations, the skewness of implied volatilities as a function of strike is jointly determined by the form of $\sigma$ and the correlation between $W_t$ and ${\alpha }_t$, whereas the convexity of implied volatilities primarily originates with the variance of ${\alpha }_t$ and, in some cases, with the convexity of $\sigma$. Popular choices for $M_t$ include the lognormal dynamics of SABR, and the mean-reverting Cox-Ingersoll-Ross (CIR) dynamics of the Heston process.

We are not focused on using (1) and (2) in a fully dynamic sense, as needed, for instance, for an arbitrage-free model for path-dependent interest rate options. Instead, we wish to use the model in the same way that the SABR model is used in practice: as a model for marking European option prices at single maturities, likely with different model parameters for each maturity. This will ultimately affect how we treat certain drift terms that arise in theory; our focus is on preserving only time 0 forwards, rather than on imposing hard intertemporal martingale conditions.

Local volatility: a reduced-form approach

In this section, we show how to reduce the pricing of the LSV model into independent components via a replication argument, and we present approximate options pricing formulas (4) and (5) in terms of one-dimensional integrals (using either call and put prices on $M_T$ or the CDF of $M_T$).

Lamperti’s harmonic transform and drift approximation.

For a positive and sufficiently regular LV function $\sigma$, we consider the Lamperti transform:

$$G(f)={\int }_{F_0}^f\frac{\mathrm{d}u}{\sigma (u)}$$

and the associated process $g_t\triangleq G(F_t)$. An application of Ito’s lemma to $g_t$ results in:

	$g_T$	$=M_T-{\mu }_T$
	${\mu }_T$	$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^T}{\alpha }_t^2{\sigma }^{\prime }(G^{-1}(g_t))\mathrm{d}t$

Intuitively, the Lamperti transform is a tool of choice to disentangle local from pure stochastic volatility: the effect of the local component is then captured via a straight change in co-ordinates, giving rise to a stochastic drift, ${\mu }_T$, that ensures the martingality of $F_T=G^{-1}(g_T)$. As explained earlier, we are solely concerned with the term distribution of $F_T$ for a given horizon $T$, so we make the pragmatic assumption that ${\mu }_T=\mu$, where $\mu$ is an ad hoc constant, fixed at a value that ensures preservation of the forward, in the sense that:

$$?(F_T)=?(G^{-1}(M_T-\mu ))=F_0$$

where $?$ is the expectation operator under the pricing measure $\mathbb{Q}$. Accordingly, we adopt the following definition for the term distribution of $F_T$ from now on:

$$F_T=G^{-1}(M_T-\mu )$$

(3)

The simplifying assumption made on the stochastic drift ${\mu }_T$ impacts only mildly the distribution of the underlying. This is certainly the case as time to maturity decreases, where $F_T$ in (3) will converge to the exact solution of the stochastic differential equation (SDE), which is consistent with Berestycki et al (2004). For longer expiries, one can expect the model to exhibit larger differences with the true LSV SDE, but qualitatively the LV parametrisation will have a similar effect on the implied volatility smile. Note that, regardless of the maturity, the model is exact for normal/lognormal dynamics and for pure SV, providing increased confidence for the general case. As per construction, (3) will produce a valid (ie, non-negative) density for $F_T$ and is not restrictive in terms of the resulting probability distribution, since the LV $\sigma$, being generic, can absorb any differences when calibrated.

Option pricing

For any given random variable $Y_T$ observed at $T$, we adopt the following notation throughout: $C^Y(K)=?{(Y_T-K)}^{+}$ (respectively, $P^Y(K)=?{(K-Y_T)}^{+}$) for a call value (respectively, a put value) and ${\mathrm{\Phi }}^Y(K)=?(1_{Y_T\le K})$ for its CDF.

Now, since $F_T=G^{-1}(M_T-\mu )$, it follows that $C^F(K)$ and $P^F(K)$ can be written in terms of call/put options on $M_T$, via a replication argument:

	$C^F(K)$	$=\sigma (K)C^M(G(K)+\mu )+{\displaystyle {\int }_K^{\mathrm{\infty }}}{\sigma }^{\prime }(k)C^M(G(k)+\mu )\mathrm{d}k$
	$P^F(K)$	$=\sigma (K)P^M(G(K)+\mu )-{\displaystyle {\int }_{-\mathrm{\infty }}^K}{\sigma }^{\prime }(k)P^M(G(k)+\mu )\mathrm{d}k$		(4)

We may also write these equations in terms of ${\varphi }^M$, the CDF of $M_T$:

	$C^F(K)$	$={\displaystyle {\int }_K^{\mathrm{\infty }}}\{1-{\varphi }^M(G(k)+\mu )\}\mathrm{d}k$
	$P^F(K)$	$={\displaystyle {\int }_{-\mathrm{\infty }}^K}{\varphi }^M(G(k)+\mu )\mathrm{d}k$		(5)

Proof (sketch).

Results in (4) are based on a replication argument using Carr-Madan integral representation. Identities in (5) stem from the classical relationship between option prices and the CDF (see Bang (2018) for details). ∎

The integrals involved in (4) and (5) can be efficiently implemented, provided we have either fast closed-form solutions for options on $M_T$ (see applications to SABR-style models) or specific forms for the LV function $\sigma$ (see applications to Heston models with piecewise linear LV).

SABR LV model

In this section, we use the normal SABR process as the SV building block $M_t$ in our model. As mentioned earlier, fast closed-form expressions for options on $M_t$ are critical to the practical implementation of (4) or (5), so we present efficient, arbitrage-free approximation formulas for options prices and the CDF for normal SABR dynamics.

Normal SABR

We recall the normal SABR model dynamics sets $F_t=M_t$, where:

	$\mathrm{d}{M}_t$	$={\alpha }_t\mathrm{d}{W}_t$
	$\mathrm{d}{\alpha }_t$	$=\nu {\alpha }_t(\rho \mathrm{d}{W}_t+\sqrt{1-{\rho }^2}\mathrm{d}{W}_t^{\perp })$

Here, $\nu >0$, $\rho \in (-1,1)$, and $W_t$ and $W_t^{\perp }$ are two independent Brownian motions. From Gyongy’s lemma, we know that there exists a unique, pure LV process ${\tilde{M}}_t$ with the same marginal distributions as $M_t$. In lemma 1, we describe explicitly the dynamics of ${\tilde{M}}_t$.

Lemma 1.

Consider the process $z_t$, defined by:

$$\mathrm{d}{z}_t=\nu \sqrt{1+2\rho z_t+z_t^2}\mathrm{d}{W}_t,z_0=0$$

Then, for any $t\mathrm{\ge }\mathrm{0}$, ${\stackrel{\mathrm{~}}{M}}_t\mathrm{=}{M}_{\mathrm{0}}\mathrm{+}\mathrm{(}{\alpha }_{\mathrm{0}}\mathrm{/}\nu \mathrm{)}\mathit{}{z}_t$ has the same marginal distributions as $M_t$.

Proof (sketch).

Given the result of Gyongy, we know the squared volatility term in the process ${\tilde{M}}_t$ is $?({\alpha }_t^2|M_t=x)$, ensuring consistent marginal distributions between ${\tilde{M}}_t$ and $M_t$. An explicit calculation of this expression can be performed by a measure change, along the lines of Balland & Tran (2013). ∎

Let us define ${\xi }_T\triangleq {\mathrm{e}}^{\nu h_T}$, where:

$$h_T={\int }_0^{z_T}\frac{\mathrm{d}z}{\nu \sqrt{1+2\rho z+z^2}}$$

Lemma 2 relates options on $M_T$ to options on ${\xi }_T$.

Lemma 2.

For a given strike $K$, let $z^K\mathrm{\triangleq }\mathrm{(}\nu \mathrm{/}{\alpha }_{\mathrm{0}}\mathrm{)}\mathit{}\mathrm{(}K\mathrm{-}{M}_{\mathrm{0}}\mathrm{)}$ and set:

$${\xi }^K\triangleq \frac{\sqrt{1-{\rho }^2+{(z^K+\rho )}^2}+z^K+\rho }{1+\rho }$$

We then have:

	$C^M(K)$	$={\displaystyle \frac{{\alpha }_0}{2\nu }}\left[(1+\rho )C^{\xi }({\xi }^K)+(1-\rho )P^{1/\xi }\left({\displaystyle \frac{1}{{\xi }^K}}\right)\right]$		(6)
	${\mathrm{\Phi }}^M(K)$	$={\mathrm{\Phi }}^{\xi }({\xi }^K)$

Proof.

Integrating explicitly the relation:

$$h_T={\int }_0^{z_T}{(\nu \sqrt{1+2\rho z+z^2})}^{-1}\mathrm{d}z$$

and inverting the resulting expression leads to the reconstruction formula:

$$z_T=-\rho +\frac{1}{2}\left((1+\rho ){\xi }_T-\frac{1-\rho }{{\xi }_T}\right)$$

(7)

From the definitions of $z^K$ and ${\xi }^K$, we have:

$$z^K=-\rho +\frac{1}{2}\left((1+\rho ){\xi }^K-\frac{1-\rho }{{\xi }^K}\right)$$

Thus, subtracting $z^K$ from $z_T$ and reverting back to the initial variable ${\tilde{M}}_T$, we see:

$$M_T-K\stackrel{d}{=}\frac{{\alpha }_0}{2\nu }\left[(1+\rho )({\xi }_T-{\xi }^K)+(1-\rho )\left(\frac{1}{{\xi }^K}-\frac{1}{{\xi }_T}\right)\right]$$

which shows $M_T-K$ is distributed as a sum of two increasing functions in ${\xi }_T$, vanishing simultaneously at ${\xi }^K$. Thus, we clearly have:

	${(M_T-K)}^{+}$	$\stackrel{d}{=}{\displaystyle \frac{{\alpha }_0}{2\nu }}\left[(1+\rho ){({\xi }_T-{\xi }^K)}^{+}+(1-\rho ){\left({\displaystyle \frac{1}{{\xi }^K}}-{\displaystyle \frac{1}{{\xi }_T}}\right)}^{+}\right]$		(8)
	$1_{M_t\le K}$	$\stackrel{d}{=}{1}_{{\xi }_t\le {\xi }^K}$		(9)

Applying the expectation operator $?$ then proves the lemma. ∎

We remark that the elegant technique used to prove (8)–(9) is known as Jamshidian’s trick, which was originally used in the context of swaption valuation in a one-factor model. Our next goal is to work out accurate and non-arbitrageable closed-form approximations for $C^{\xi }$, $P^{1/\xi }$ and ${\mathrm{\Phi }}^{\xi }$. We start by using Ito’s lemma to note that:

$$\mathrm{d}{h}_t=\mathrm{d}{W}_t-\frac{\nu }{2}\tanh (\nu (h_t+\overline{h}))\mathrm{d}t$$

(10)

with:

$$\overline{h}\triangleq \frac{1}{2\nu }\ln \left(\frac{1+\rho }{1-\rho }\right)$$

The dynamics of $h_t$ in (10) resemble a mean-reverting drifted Brownian motion, so both its mean and variance will be affected by the (locally) stochastic drift term. Given ${\xi }_T$ is defined as an exponential of $h_t$, one is tempted to use a simple lognormal approximation for ${\xi }_T$, with moment-matched mean and variance. Indeed, we find that this works quite well but nonetheless present an alternative in proposition 1, which boasts even better precision. The idea behind the proposition is to apply a measure change to capture the subtle dynamics of $h_T$ more accurately, without sacrificing tractability.

Proposition 1.

The following (arbitrage-free) expressions can be used as very accurate closed-form solutions for option prices and the CDF for the process ${\xi }_T$:

	$C^{\xi }({\xi }^K)$	$\approx {\displaystyle \frac{3\mathrm{\Gamma }}{2\gamma }}C(\nu ,\nu ,k,\overline{h},T)-{\displaystyle \frac{\mathrm{\Gamma }}{2\gamma }}C(\nu ,3\nu ,k,\overline{h},T)$		(11)
	$P^{1/\xi }\left({\displaystyle \frac{1}{{\xi }^K}}\right)$	$\approx {\displaystyle \frac{1}{2\gamma \mathrm{\Gamma }}}C(-\nu ,3\nu ,k,\overline{h},T)-{\displaystyle \frac{3}{2\gamma \mathrm{\Gamma }}}C(-\nu ,\nu ,k,\overline{h},T)$		(12)
	$1-{\mathrm{\Phi }}^{\xi }({\xi }^K)$	$\approx {\displaystyle \frac{3}{2\gamma }}\mathrm{{\rm Y}}(\nu ,k,\overline{h},T)-{\displaystyle \frac{1}{2\gamma }}\mathrm{{\rm Y}}(3\nu ,k,\overline{h},T)$		(13)

with:

	$k$	$={\displaystyle \frac{1}{\nu }}\ln \left({\displaystyle \frac{{\xi }^K}{\mathrm{\Gamma }}}\right)$		(14)
	$\gamma$	$=\frac{3}{2}{\mathrm{{\rm Y}}}^{\star }(\nu ,\overline{h},T)-\frac{1}{2}{\mathrm{{\rm Y}}}^{\star }(3\nu ,\overline{h},T)$		(15)
	$\mathrm{\Gamma }$	$={\displaystyle \frac{\rho +\sqrt{{\rho }^2+(1-{\rho }^2)E^{+}{E}^{-}}}{(1+\rho )E^{+}}}$		(16)
	$E^{\pm }$	$={\displaystyle \frac{{\mathrm{e}}^{{\nu }^{2}T/2}}{2\gamma }}(3{\mathrm{{\rm Y}}}^{\star }(\nu ,\overline{h}\pm \nu T,T)-{\mathrm{{\rm Y}}}^{\star }(3\nu ,\overline{h}\pm \nu T,T))$		(17)

and where the auxiliary functions ${\mathrm{{\rm Y}}}^{\mathrm{\star }}$, $\mathrm{{\rm Y}}$ and $C$ have been defined below via the Gaussian kernel and re-expressed, using elementary integral manipulations, in terms of the normal CDF, $N\mathit{}\mathrm{(}\mathrm{\cdot }\mathrm{)}$:

	${\mathrm{{\rm Y}}}^{\star }(\lambda ,\overline{h},T)$	$\triangleq {\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}{\displaystyle \frac{{\mathrm{e}}^{-(\lambda /2)|h+\overline{h}|}}{\sqrt{2\pi T}}}{\mathrm{e}}^{-h^2/2T}\mathrm{d}h$
		$={\mathrm{e}}^{{\lambda }^{2}T/8}\left[{\mathrm{e}}^{-\lambda \overline{h}/2}N\left({\displaystyle \frac{\overline{h}}{\sqrt{T}}}-{\displaystyle \frac{\lambda \sqrt{T}}{2}}\right)+{\mathrm{e}}^{\lambda \overline{h}/2}N\left(-{\displaystyle \frac{\overline{h}}{\sqrt{T}}}-{\displaystyle \frac{\lambda \sqrt{T}}{2}}\right)\right]$
	$C(\sigma ,\lambda ,k,\overline{h},T)$	$\triangleq {\displaystyle {\int }_k^{+\mathrm{\infty }}}{({\mathrm{e}}^{\sigma h}-{\mathrm{e}}^{\sigma k})}^{+}{\displaystyle \frac{{\mathrm{e}}^{-(\lambda /2)|h+\overline{h}|}}{\sqrt{2\pi T}}}{\mathrm{e}}^{-h^2/2T}\mathrm{d}h$
		$={\mathrm{e}}^{{\sigma }^{2}T/2}\mathrm{{\rm Y}}(\lambda ,k-\sigma T,\overline{h}+\sigma T,T)-{\mathrm{e}}^{\sigma k}\mathrm{{\rm Y}}(\lambda ,k,\overline{h},T)$		(18)

with:

	$\mathrm{{\rm Y}}(\lambda ,k,\overline{h},T)$	$\triangleq {\displaystyle {\int }_k^{+\mathrm{\infty }}}{\displaystyle \frac{{\mathrm{e}}^{-(\lambda /2)|h+\overline{h}|}}{\sqrt{2\pi T}}}{\mathrm{e}}^{-h^2/2T}\mathrm{d}h$
		$={\mathrm{e}}^{({\lambda }^{2}T/8)-(\lambda \overline{h}/2)}N\left(-{\displaystyle \frac{k}{\sqrt{T}}}-{\displaystyle \frac{\lambda \sqrt{T}}{2}}\right)\mathit{\hspace{1em}}\mathit{\text{if }}k\ge -\overline{h}$
				(19)

Proof.

We consider the exponential martingale ${\theta }_t^{\star }$ in the measure $\mathbb{Q}$:

$$\frac{\mathrm{d}{\theta }_t^{\star }}{{\theta }_t^{\star }}=\frac{\nu }{2}\tanh (\nu (h_t+\overline{h}))\mathrm{d}{W}_t$$

The associated measure ${\mathbb{Q}}^{\star }$ is defined via:

$$\frac{\mathrm{d}{\mathbb{Q}}^{\star }}{\mathrm{d}\mathbb{Q}}={\theta }_T^{\star }$$

By construction, $h_T$ is a driftless Brownian motion under $Q^{\star }$. Denoting by ${\delta }_h$ the Dirac distribution at $h$:

$$?({\delta }_h(h_T))={?}^{\star }\left(\frac{{\delta }_h(h_T)}{{\theta }_T^{\star }}\right)={?}^{\star }\left({?}^{\star }[\frac{1}{{\theta }_T^{\star }}|h_T]{\delta }_h(h_T)\right)$$

where ${?}^{\star }$ denotes expectation in measure ${\mathbb{Q}}^{\star }$. Defining:

$$f(h)=\sqrt{\cosh (\nu (h+\overline{h}))}$$

and applying Ito’s lemma, we find $\mathrm{d}(f(h_t)/{\theta }_t^{\star })=O(\mathrm{d}t)$. Moreover, using a second-order Taylor expansion in ${\mathrm{e}}^{-\nu |h+\overline{h}|}$, we see:

$$\frac{1}{f(h)}\approx {\mathrm{e}}^{-(\nu /2)|h+\overline{h}|}\frac{3-{\mathrm{e}}^{-\nu |h+\overline{h}|}}{2}$$

Thus, we can write:

	${?}^{\star }[{\displaystyle \frac{1}{{\theta }_T^{\star }}}|h_T]$	$={\displaystyle \frac{1}{f(h_T)}}{?}^{\star }[{\displaystyle \frac{f(h_T)}{{\theta }_T^{\star }}}|h_T]$
		$\approx {\displaystyle \frac{1}{\gamma }}{\mathrm{e}}^{-(\nu /2)|h_T+\overline{h}|}{\displaystyle \frac{3-{\mathrm{e}}^{-\nu |h_T+\overline{h}|}}{2}}\triangleq {\displaystyle \frac{1}{{\theta }^{†}(h_T)}}$

The constant $\gamma$ is introduced to ensure:

$${?}^{Q^{\star }}\left(\frac{1}{{\theta }^{†}(h_T)}\right)=1$$

Finally, we can write ${\xi }_T={\mathrm{e}}^{\nu h_T}\approx \mathrm{\Gamma }{\mathrm{e}}^{\nu h_T^{†}}$, where the density of $h_T^{†}$ under $Q$ is:

$$\varrho (h)=\frac{{\mathrm{e}}^{-(\nu /2)|h+\overline{h}|}(3-{\mathrm{e}}^{-\nu |h+\overline{h}|})}{\mathrm{2\gamma }}\frac{{\mathrm{e}}^{-h^2/2T}}{\sqrt{2\pi T}}$$

(20)

$\gamma$ can be found by enforcing ${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\varrho (h)\mathrm{d}h=1$, which leads, after standard integral manipulations, to the expression in (15).

The parameter $\mathrm{\Gamma }$, which has been introduced to account for the small measure change approximations made on the way, can be found using the martingale condition $?(z_T)=0$ in (7), that is:

$$-\rho +\frac{1}{2}\left((1+\rho )\mathrm{\Gamma }{E}^{+}-\frac{1-\rho }{\mathrm{\Gamma }}{E}^{-}\right)=0$$

where:

$$E^{\pm }\triangleq {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\mathrm{e}}^{\pm \nu h}\varrho (h)\mathrm{d}h$$

can be conveniently expressed as in (17). Thus, $\mathrm{\Gamma }$ is the solution of a second-order equation with a unique admissible solution:

$$\mathrm{\Gamma }=\frac{\rho +\sqrt{{\rho }^2+(1-{\rho }^2)E^{+}{E}^{-}}}{(1+\rho )E^{+}}$$

Note (20) is a valid density by construction and completely tractable. As a matter of fact, since ${\xi }_T=\mathrm{\Gamma }{\mathrm{e}}^{\nu h_T}$, and from the definition of $k$ in (14), we have:

$C^{\xi }({\xi }^K)$

$\triangleq ?{({\xi }_T-{\xi }^K)}^{+}\approx \mathrm{\Gamma }{\displaystyle {\int }_k^{+\mathrm{\infty }}}({\mathrm{e}}^{\nu h}-{\mathrm{e}}^{\nu k})\varrho (h)\mathrm{d}h$

which can be expressed as in (11) using the function $C$ defined in (18). The proof for (12) is similar. For (13), we note that:

$$1-{\mathrm{\Phi }}^{\xi }({\xi }^K)=?(1_{{\xi }_T>{\xi }^K})=?(1_{h_T>k})\approx {\int }_k^{+\mathrm{\infty }}\varrho (h)\mathrm{d}h$$

which yields the result. ∎

Bringing together results from lemma 2 and proposition 1, we arrive at explicit and fast approximations for options on $M_T$. We emphasise that these expressions are arbitrage-free by construction, for any set of parameters and for any maturity. Empirically, the approximations perform surprisingly well, even for large maturities, as is shown by the results in figure 1; this compares our approximation to both the Hagan expansion in Hagan et al (2002) and to exact values in Antonov et al (2017) (computed by an expensive integral-based closed-form solution). The new formulas only require a few calls to the Gaussian CDF (with a cost similar to that of the Hagan formula).

Figure 1: Comparison of implied probability density function (left $y$-axis) and normal implied volatility (right $y$-axis) for our new approach (‘new’), the exact approach in Antonov et al (2017) (‘exact’) and the asymptotic approximation in Hagan et al (2002) (‘Hagan’) for the normal SABR model. The swaption tenor is 10 years for a 50-year maturity. Forward $F_{\text{0}}=\text{2}\%$, at-the-money normal volatility $\widehat{\sigma }=\text{0.5}\%$, volatility of volatility $\nu =\text{20}\%$ and correlation $\rho =\text{50}\%$.

Local volatility specification

Now the problem of determining fast and accurate approximations for the normal SABR model is out of the way, we can freely design the LV function $\sigma$ in (1) to provide trading desks, using a model that can be fine-tuned to accommodate their specific requirements, for instance:

• •

  adequate support, eg, the distribution for $F_T$ should allow for negative rates;

• •

  flexibility, the parametrisation needs to be generic enough to be calibrated to a reasonably large number of liquid option strikes (eg, five);

• •

  control over the dynamics between level and volatility of the underlying for risk computation purposes (backbone); this should comply with trader requirements and empirical evidence; and

• •

  control over model behaviour at high and low strikes and, ultimately, over convex products such as CMS.

While many different LV specifications could be used to target these requirements, let us give a reasonable example for practical purposes. To address the first point, we choose an LV function that remains positive for negative values of $F$. To this aim, for $\epsilon \ge 0$, we introduce the function $\eta (\epsilon ,f)=\epsilon \ln (1+{\mathrm{e}}^{f/\epsilon })$, which regularises $f^{+}$. To control the backbone, we consider a constant elasticity of variance (power) specification, with level-dependent power:

$$\beta (f)=\frac{1}{2}({\beta }_l+{\beta }_h)+\frac{1}{2}({\beta }_h-{\beta }_l)\tanh (\lambda (f-F_0))$$

Note $\beta (f)$ varies smoothly from ${\beta }_l$ to ${\beta }_h$, with a speed controlled by the constant $\lambda$. So equipped, we are now able to explicitly define the LV function:

$${\sigma }^{\star }(f)={(\eta (\epsilon ,f))}^{\beta (f)}+\psi \eta ({\epsilon }_s,f-K_h)$$

In the function ${\sigma }^{\star }$, the parameter $\epsilon$ controls the smile for low strikes and typically the probability mass below 0. If $\epsilon =0$, then $F_T$ remains above 0 (if it started above 0). $K_h$ is an arbitrary strike outside the region of observable strikes, but beyond which the smile still affects CMS forwards. The parameter $\psi$ controls the volatility beyond $K_h$, and the parameter ${\epsilon }_s$ is used as a smoothing parameter. To avoid explosive behaviour from ${\sigma }^{\star }$ in the wings, we finally choose to extrapolate the LV flat outside some interval $[\underset{¯}{f},\overline{f}]$. One benefit of this extrapolation is the integrals in (4) are truncated over finite support, making them more amenable to fast numerical integration. For interest rate applications, a possible parameterisation for $\underset{¯}{f}$, $\overline{f}$ and $K_h$ is shown in figure 2. The final expression for the LV reads:

$\sigma (f)={\sigma }^{\star }(\min (\max (f,\underset{¯}{f}),\overline{f}))$

(21)

In figure 1 we demonstrate the impact of the parameters controlling the wings of the LV.

Figure 2: Control over the wings of the LV (normalised at the forward $F_{\text{0}}=\text{0.02}$), varying $\epsilon$ and $\psi$. Base case (square-root LV) corresponds to ${\beta }_l={\beta }_h=\frac{\text{1}}{\text{2}}$, $\epsilon =\text{0}$, ${\epsilon }_s=\text{0.01}$, $\underset{¯}{f}=-\text{3}\%$, $\overline{f}=F_{\text{0}}+\text{50}\%$, $K_h=F_{\text{0}}+\text{8}\%$ and $\psi =\text{0}$.

Model settings and numerical results

We consider the LV function $\sigma$ defined in (21) via a set of LV parameters and the associated Lamperti transform $G(f)={\int }_{F_0}^f\sigma {(u)}^{-1}\mathrm{d}u$. $G(f)$ can be computed and cached, once and for all, using, for instance, a trapezoid quadrature rule on a fine grid. Options on $F_T$ struck at $K$ are then expressed in (4), where $C^M$ is computed using lemma 2 and proposition 1; put options $P^M$ may be computed using put-call parity.

For some fixed maturity $T$, our first task is calibrating $\mu$ and ${\alpha }_0$ to simultaneously reprice the forward (or, equivalently, ensuring $C^F(F_0)=P^F(F_0)$) and the at-the-money (ATM) option, $C^F(F_0)=\widehat{\sigma }\sqrt{T/(2\pi )}$, where $\widehat{\sigma }$ is the normal (Bachelier) ATM implied volatility. This is achieved using ATM call and put values, derived from setting $K=F_0$ in (4). Note that, once calibrated, the method allows us to price a strip of options simultaneously by adequately discretising the integrals in (4). In the next paragraph, the model is (successfully) calibrated to market quotes (a strip of options and a CMS forward).

Calibration to market instruments

In this section, our model is used with a full set of LV and SV parameters, calibrated to best-fit option prices and CMS forwards seen in market quotes. In figure 3, we display the probability density function (PDF) and the normal (Bachelier) volatility implied by the model, using the normal SABR formula (6) and the full LV function (21). Evidently, a very accurate calibration to these instruments can be achieved, with the resulting implied PDF being smooth and well behaved.

Figure 3: Implied PDF (left $y$-axis), LV and implied normal volatility (right $y$-axis) for a model calibrated to swaptions implied volatilities (Market), and the CMS convexity adjustments (61 basis points). Forward $F_{\text{0}}=\text{1.833}\%$, ATM volatility $\widehat{\sigma }=\text{0.533}\%$, volatility of volatility $\nu =\text{16.5}\%$, correlation $\rho =\text{7}\%$, ${\beta }_l=\text{26.5}\%$, ${\beta }_h=\text{0}\%$, $\epsilon =\text{0.0032}$, ${\epsilon }_s=\text{0.001}$, $\lambda =\text{0.5}$ and $\psi =\text{47.7}$.

Heston model with linear piecewise local volatility

As mentioned earlier, the proposed method for building LSV term volatility models is generic and can accommodate different dynamics. In this section, we briefly illustrate the breadth of the model, using it to equip the Heston model with a general piecewise linear LV. For this purpose, we consider, as a building block, a normal Heston dynamics for the process $M_t$:

	$\mathrm{d}{M}_t$	$=\sqrt{{\zeta }_t}\mathrm{d}{W}_t$
	$\mathrm{d}{\zeta }_t$	$=\kappa (\theta -{\zeta }_t)\mathrm{d}t+\nu \sqrt{{\zeta }_t}(\rho \mathrm{d}{W}_t+\sqrt{1-{\rho }^2}\mathrm{d}{W}^{\perp })$

where ${\zeta }_0$, $\rho$ and $\nu$ control the ATM level, the skew and the convexity of the implied volatility of the model, respectively.

The characteristic function of $M_T$, $\chi (u)=E({\mathrm{e}}^{\mathrm{i}u{M}_T})$, is known in closed form; thus, the Lewis-Lipton formula can be employed to recover option prices on $M_T$ or on ${\mathrm{e}}^{\lambda M_T}$ via an inverse Fourier integral (see Lipton 2002). The general pricing equation (4) would not be computationally workable, as it would involve double integrals. Nevertheless, in the next section we show how the assumption of piecewise linearity on the LV allows us to write options prices as a discrete sum of options in the standard Heston model.

For later use, let us denote:

$$V_{\mathrm{c},\mathrm{p}}^H(s,G)=?{\left({\epsilon }_{\mathrm{c},\mathrm{p}}\frac{e^{s(M_T-\mu -G)}-1}{s}\right)}^{+}$$

with the convention $V_{\mathrm{c},\mathrm{p}}^H(0,G)=E{({\epsilon }_{\mathrm{c},\mathrm{p}}(M_T-\mu -G))}^{+}$, ${\epsilon }_{\mathrm{c}}=1$ (respectively, ${\epsilon }_{\mathrm{p}}=-1$) for a call (respectively, for a put).

Replication formula using calls and puts

Now, assuming a piecewise linear LV (extrapolated flat, for concreteness), we consider a strip of $n$ strikes $K_{1\le i\le n}$; without loss of generality, we can assume the forward $F_0=K_{\mathrm{p}}$, with $p\in ]1,n[$. We also define the associated volatilities ${\sigma }_{1\le i\le n}$. Let:

$${\sigma }_i^{\prime }=\frac{{\sigma }_{i+1}-\sigma i}{K_{i+1}-Ki}$$

be the LV slopes, and then define the LV function:

We also denote $G_i=G(K_i)$ as the (Lamperti) transformed variables on the strike grid. The following identities hold:

	$P^F(K_1)$	$={\sigma }_1{V}_{\mathrm{p}}^H(0,G_1)$
	$P^F(K_{i+1})-P^F(K_i)$	$={\sigma }_{i+1}{V}_{\mathrm{p}}^H({\sigma }_i^{\prime },G_{i+1})-{\sigma }_i{V}_{\mathrm{p}}^H({\sigma }_i^{\prime },G_i),$
	$C^F(K_i)-C^F(K_{i+1})$	$={\sigma }_i{V}_{\mathrm{c}}^H({\sigma }_i^{\prime },G_i)-{\sigma }_{i+1}{V}_{\mathrm{c}}^H({\sigma }_i^{\prime },G_{i+1})$
	$C^F(K_n)$			(22)

Proof (sketch).

Identities in (22) can be proven by first writing the Lamperti transform $G$ (and its inverse) on each subinterval where the local volatility $\sigma$ is linear, before applying the expectation operator. Details can be found in Bang (2018). ∎

Calibration

The calibration problem involves the constraints of repricing the forward and reproducing the ATM volatility $\widehat{\sigma }$. The two parameters available for this task are the drift $\mu$ and the initial value of the volatility, ${\zeta }_0$:

	$P^F(F_0)$	$={\sigma }_1{V}_{\mathrm{p}}^H(0,G_1)+{\displaystyle \sum _{i=1}^{p-1}}{\sigma }_{i+1}{V}_{\mathrm{p}}^H({\sigma }_i^{\prime },G_{i+1})-{\sigma }_i{V}_{\mathrm{p}}^H({\sigma }_i^{\prime },G_i)$
	$C^F(F_0)$	$={\sigma }_n{V}_{\mathrm{c}}^H(0,G_n)+{\displaystyle \sum _{i=p}^{n-1}}{\sigma }_i{V}_{\mathrm{c}}^H({\sigma }_i^{\prime },G_i)-{\sigma }_{i+1}{V}_{\mathrm{c}}^H({\sigma }_i^{\prime },G_{i+1})$

Thus, the calibration problem can be written as the system:

	$P^F(F_0)-C^F(F_0)$	$=0$
	$P^F(F_0)+C^F(F_0)$	$=\widehat{\sigma }\sqrt{{\displaystyle \frac{2T}{\pi }}}$

which can be solved in $\mu$ and ${\zeta }_0$ by any standard two-dimensional root search algorithm.

Option pricing

Consider a call struck at $K\in [K_i,K_{i+1}]$. We have:

$$C^F(K)=C^F(K_{i+1})+\sigma (K)V_{\mathrm{c}}^H({\sigma }_i^{\prime },G(K))-{\sigma }_{i+1}{V}_{\mathrm{c}}^H({\sigma }_i^{\prime },G_{i+1})$$

Note $C^F(K_{i+1})$ has already been computed at the calibration stage. Thus, it is straightforward to price an option on $F_T$ struck at $K$. Note that in the degenerate case where the normal Heston process $M_t$ is Gaussian ($\nu =0$), the model can be used as an arbitrage-free LV interpolator (see Schlenkrich (2018) for details).

Conclusion

We have presented a new method for term distribution modelling that allows the user to mix an arbitrary LV with an arbitrary SV in a tractable, intuitive and arbitrage-free way. We have applied the method to the SABR models family, pairing a flexible LV with a normal SABR SV model. In the process, we developed new, efficient closed-form option pricing formulas for the normal SABR SV model. The resulting LSV model is smooth and arbitrage free, and it successfully calibrates to market quotes for swaptions and CMS forwards. Finally, we have also shown how the method can be used in the context of the Heston model by providing closed-form pricing solutions when the LV is a general linear piecewise function and the SV is a normal Heston model.

Dominique Bang is head of interest rate vanilla analytics at Bank of America Merrill Lynch in London. He is indebted to Leif Andersen and Guillaume Blacher for insightful suggestions, comments and support, and he is very grateful to Mark Lake for an efficient implementation of the model. He also thanks Guillaume Royer, Julien Pantz, Elias Daboussi and an anonymous referee for constructive feedback.

Email: [e-mail].

References

 

• Andersen L and V Piterbarg, 2010
  Interest Rate Modeling, Volume 2
  Atlantic Financial Press
• Andreasen J and B Huge, 2013
  Expanded forward volatility
  Risk January, 101–107
• Antonov A, M Konikov and M Spector, 2017
  Mixing SABR models for negative rates
  Risk September, 86–91
• Balland P and Q Tran, 2013
  SABR goes normal
  Risk May, 76–81
• Bang D, 2018
  Local-stochastic volatility models for vanilla modeling: a tractable and arbitrage free approach to option pricing
  SSRN preprint, available at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3171877
• Berestycki H, J Busca and I Florent, 2004
  Computing the implied volatility in stochastic volatility models
  Communications in Pure and Applied Mathematics 57, pages 1352–1373
• Blacher G, 2001
  A new approach for designing and calibrating stochastic volatility models for optimal delta vega hedging of exotic options
  Working Paper, Global Derivatives Conference, Juan-les-Pins, France
• Hagan P, D Kumar, A Lesniewski and D Woodward, 2002
  Managing smile risk
  Wilmott January, 84–108
• Lipton A, 2002
  The vol smile problem
  Risk February, 81–85
• Schlenkrich S, 2018
  Approximate local volatility model for vanilla rates options
  SSRN preprint, available at https://ssrn.com/abstract=3150689
• Ticot Y and X Charvet, 2013
  Inflation derivatives – LPI swaps with a smile
  Risk June, 66–71