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

## 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 ${\stackrel{~}{M}}_{t}$ with the same marginal distributions as ${M}_{t}$. In lemma 1, we describe explicitly the dynamics of ${\stackrel{~}{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 ${\stackrel{~}{M}}_{t}$ is $?({\alpha}_{t}^{2}|{M}_{t}=x)$, ensuring consistent marginal distributions between ${\stackrel{~}{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 ${\stackrel{~}{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}\mathrm{tanh}(\nu ({h}_{t}+\overline{h}))\mathrm{d}t$$ | (10) |

with:

$$\overline{h}\triangleq \frac{1}{2\nu}\mathrm{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}}\mathrm{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}\mathrm{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{\mathrm{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}^{\u2020}({h}_{T})}}$ |

The constant $\gamma $ is introduced to ensure:

$${?}^{{Q}^{\star}}\left(\frac{1}{{\theta}^{\u2020}({h}_{T})}\right)=1$$ |

Finally, we can write ${\xi}_{T}={\mathrm{e}}^{\nu {h}_{T}}\approx \mathrm{\Gamma}{\mathrm{e}}^{\nu {h}_{T}^{\u2020}}$, where the density of ${h}_{T}^{\u2020}$ 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).

### 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 \mathrm{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})\mathrm{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{\xaf}{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{\xaf}{f}$, $\overline{f}$ and ${K}_{h}$ is shown in figure 2. The final expression for the LV reads:

$\sigma (f)={\sigma}^{\star}(\mathrm{min}(\mathrm{max}(f,\underset{\xaf}{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{\xaf}{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) |

### 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: dominique.bang@baml.com.

## 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

Only users who have a paid subscription or are part of a corporate subscription are able to print or copy content.

To access these options, along with all other subscription benefits, please contact info@risk.net or view our subscription options here: http://subscriptions.risk.net/subscribe

You are currently unable to print this content. Please contact info@risk.net to find out more.

You are currently unable to copy this content. Please contact info@risk.net to find out more.

Copyright Infopro Digital Limited. All rights reserved.

You may share this content using our article tools. Printing this content is for the sole use of the Authorised User (named subscriber), as outlined in our terms and conditions - https://www.infopro-insight.com/terms-conditions/insight-subscriptions/

If you would like to purchase additional rights please email info@risk.net

Copyright Infopro Digital Limited. All rights reserved.

You may share this content using our article tools. Copying this content is for the sole use of the Authorised User (named subscriber), as outlined in our terms and conditions - https://www.infopro-insight.com/terms-conditions/insight-subscriptions/

If you would like to purchase additional rights please email info@risk.net

####
More on Derivatives

##### One strike and they’re out: traders threaten liquidity stoppage

PTFs vow to withdraw from Treasuries market in protest at SEC registration plans

##### Collateral markets in need of rewiring

New data suggests a tech upgrade is needed to avoid a large central bank footprint in markets

##### Dutch pensions have extra year to restructure hedges

January 2028 implementation date allows more time for long-dated swaps to roll off

##### Fast LPs accuse rivals of maxing out last look response times

Firms with sub-10ms checks complain of losing volumes to slower rivals, prompting one to ditch ECNs

##### Swap Connect shines light on US client clearing hurdles

New scheme may intensify calls for CFTC to reassess its exempt DCO limitations

##### US life insurer index options market hits $1trn mark

Counterparty Radar: Lincoln Financial emerges as top player in Q4 with $43 billion portfolio increase

##### Long-end euro swap pricing anomaly remains largely untapped

Deviation in swap curve attracts limited interest because of regulatory and pension reform barriers

##### SG1 growth slower than expected, say LPs

Despite sluggish take-up of Singapore FX matching engines, some hope a new NDF venue will offer a boost