Fast calibration of the Libor Market Model with Stochastic Volatility and Displaced Diffusion

TL;DR

This paper introduces a fast calibration method for the Libor Market Model with Stochastic Volatility and Displaced Diffusion using Edgeworth and Gram-Charlier expansions, significantly reducing computational time.

Contribution

It combines analytical moment-based density approximations with existing models to enhance calibration efficiency and interpretability.

Findings

Achieves 98% reduction in calibration time

Provides a smile formula relating volatility to moneyness

Demonstrates effectiveness of Edgeworth and Gram-Charlier expansions

Abstract

This paper demonstrates the efficiency of using Edgeworth and Gram-Charlier expansions in the calibration of the Libor Market Model with Stochastic Volatility and Displaced Diffusion (DD-SV-LMM). Our approach brings together two research areas; first, the results regarding the SV-LMM since the work of Wu and Zhang (2006), especially on the moment generating function, and second the approximation of density distributions based on Edgeworth or Gram-Charlier expansions. By exploring the analytical tractability of moments up to fourth order, we are able to perform an adjustment of the reference Bachelier model with normal volatilities for skewness and kurtosis, and as a by-product to derive a smile formula relating the volatility to the moneyness with interpretable parameters. As a main conclusion, our numerical results show a 98% reduction in computational time for the DD-SV-LMM calibration process compared to the classical numerical integration method developed by Heston (1993).

Tables

\@tabular@row@before@xcolor   \@xcolor@tabular@before Method	CPU Time (seconds)
\@tabular@row@before@xcolor   \@xcolor@row@afterHeston	425.1
\@tabular@row@before@xcolor   \@xcolor@row@afterEdgeworth	8.2
\@tabular@row@before@xcolor   \@xcolor@row@after

Equations

$$F_j(t)=\frac{1}{\Delta T_j}\left(\frac{P(t,T_j)}{P(t,T_{j+1})}-1\right).$$

$$\text{ for }t\le T_j,\mathrm{d}{F}_j(t)=(F_j(t)+\delta ){\zeta }_j(t)\cdot \mathrm{d}{Z}_t^{j+1},$$

$$\mathrm{d}V(t)=\kappa \left(\theta -V(t)\right)\mathrm{d}t+ϵ\sqrt{V(t)}\mathrm{d}{W}_t,$$

$$\mathrm{d}{F}_j(t)=(F_j(t)+\delta )\sqrt{V(t)}{\gamma }_j(t)\cdot \left(\mathrm{d}{Z}_t-{\sigma }_{j+1}(t)\sqrt{V(t)}\mathrm{d}t\right),$$

$${\sigma }_{j+1}(t)=-\sum _{k=m(t)}^j\frac{\Delta T_k(F_k(t)+\delta )}{1+\Delta T_k{F}_k(t)}{\gamma }_k(t),$$

$${\rho }_j(t)\mathrm{d}t=\mathbb{E}\left[\left(\frac{{\gamma }_j(t)}{\Vert {\gamma }_j(t)\Vert }\cdot \mathrm{d}{Z}_t\right)\mathrm{d}{W}_t\right].$$

$$R_{m,n}(t)=\frac{P(t,T_m)-P(t,T_n)}{B^S(t)},$$

$$PS(0,K)=B^S(0){\mathbb{E}}^S\left[\max (R_{m,n}(T_m)-K,0)\right].$$

$$\begin{array}{rl}\mathrm{d}{R}_{m,n}(t)& =\sqrt{V(t)}\sum _{j=m}^{n-1}\frac{\partial R_{m,n}(t)}{\partial F_j}(F_j(t)+\delta ){\gamma }_j(t)\cdot \mathrm{d}{Z}_t^S,\\ \mathrm{d}V(t)& =\kappa \left(\theta -{\tilde{\xi }}^S(t)V(t)\right)\mathrm{d}t+ϵ\sqrt{V(t)}\mathrm{d}{W}_t^S,\end{array}$$

$$\frac{\partial R_{m,n}(t)}{\partial F_j}={\alpha }_j(t)+\frac{\Delta T_j}{1+\Delta T_j{F}_j(t)}\sum _{k=m}^{j-1}{\alpha }_k(t)\left(F_k(t)-R_{m,n}(t)\right).$$

$$w_j(0)=\frac{\partial R_{m,n}(0)}{\partial F_j}(F_j(0)+\delta ),$$

$$\begin{array}{rl}\mathrm{d}{R}_{m,n}(t)& =\sqrt{V(t)}\sum _{j=m}^{n-1}{w}_j(0){\gamma }_j(t)\cdot \mathrm{d}{Z}_t^S,0\le t<T_m,\\ \mathrm{d}V(t)& =\kappa \left(\theta -{\tilde{\xi }}_0^S(t)V(t)\right)\mathrm{d}t+ϵ\sqrt{V(t)}\mathrm{d}{W}_t^S,\end{array}$$

$$\psi \left(R_{m,n}(t),V(t),t;z\right)={\mathbb{E}}^S\left[e^{z{R}_{m,n}(T_m)}|{\mathcal{F}}_t\right],\text{leavevmode}\text{leavevmode}\text{leavevmode}z\in \mathbb{R}.$$

$$\frac{\partial \psi }{\partial t}+\left(\kappa \theta -\kappa \xi V\right)\frac{\partial \psi }{\partial V}+\frac{1}{2}{ϵ}^{2}V\frac{{\partial }^2\psi }{{\partial V}^2}+ϵ\rho \lambda V\frac{{\partial }^2\psi }{\partial V\partial x}+\frac{1}{2}{\lambda }^{2}V\frac{{\partial }^2\psi }{{\partial x}^2}=0,$$

$$\xi \equiv {\tilde{\xi }}_0^S(t),\lambda \equiv \Vert \sum _{j=m}^{n-1}{w}_j(0){\gamma }_j(t)\Vert \text{ and }\rho =\frac{1}{\lambda }\sum _{j=m}^{n-1}{w}_j(0)\Vert {\gamma }_j(t)\Vert {\rho }_j(t),$$

$$\psi \left(x,V,t;z\right)=e^{A(\tau ,z)+B(\tau ,z)V+zx},$$

$$\{\begin{array}{rl}\frac{\partial A}{\partial \tau }& =\kappa \theta B,\\ \frac{\partial B}{\partial \tau }& =\frac{1}{2}{ϵ}^2{B}^2+\left(ϵ\rho \lambda z-\kappa \xi \right)B+\frac{1}{2}{\lambda }^2{z}^2,\end{array}$$

$$\{\begin{array}{rl}A(\tau ,z)& =A\left({\tau }_j,z\right)+{\tilde{A}}_j(\tau ,z)\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall \tau \in ({\tau }_j,{\tau }_{j+1}],\\ B(\tau ,z)& =B\left({\tau }_j,z\right)+{\tilde{B}}_j(\tau ,z)\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall \tau \in ({\tau }_j,{\tau }_{j+1}],\end{array}$$

$$f(z)=\phi (z)\sum _{n=0}^{\infty }{c}_n{H}_n(z),$$

$$H_n(z)\phi (z)={\phi }^{(n)}(z).$$

$$Z=\frac{X-\mathbb{E}[X]}{\nu }.$$

$$g_1(z)=\phi (z)\left\{1-\frac{{\mu }_3}{6}{H}_3(z)+\frac{{\mu }_4-3}{24}{H}_4(z)\right\}.$$

$$\mathbb{P}(S_n\le x)\approx \Phi (x)-\frac{{\gamma }_3}{6\sqrt{n}}\phi (x)H_2(x)+\frac{{\gamma }_4-3}{24n}\phi (x)H_3(x)+\frac{{\gamma }_3^2}{72n}\phi (x)H_5(x).$$

$$\mathbb{P}(Z\le z)\approx \Phi (z)-\frac{{\mu }_3}{6}\phi (z)H_2(z)+\frac{{\mu }_4-3}{24}\phi (z)H_3(z)+\frac{{\mu }_3^2}{72}\phi (z)H_5(z).$$

$$g_2(z)=g_1(z)+\phi (z)\frac{{\mu }_3^2}{72}{H}_6(z),$$

$$Z=\frac{R_{m,n}(T_m)-R_{m,n}(0)}{\nu },$$

$$PS(0,K)=B^S(0){\mathbb{E}}^S\left[\max (R_{m,n}(0)+\nu Z-K,0)\right].$$

$$\widetilde{PS}(0,K)=B^S(0){\int }_{z_K}^{\infty }\left(R_{m,n}(0)+\nu z-K\right)g(z)\mathrm{d}z=\nu B^S(0){\int }_{z_K}^{\infty }\left(z-z_K\right)g(z)\mathrm{d}z.$$

$${\widetilde{PS}}_0(0,K)=\nu B^S(0)\left\{\phi (z_K)-z_K\Phi (-z_K)\right\},$$

$${\widetilde{PS}}_1(0,K)={\widetilde{PS}}_0(0,K)+\nu B^S(0)\phi (z_K)\left\{\frac{{\mu }_3}{6}{z}_K+\frac{{\mu }_4-3}{24}(z_k^2-1)\right\},$$

Full text

11footnotetext: Email: laurent.devineau@milliman.com22footnotetext: Email: pierre-edouard.arrouy@milliman.com33footnotetext: Email: paul.bonnefoy@milliman.com44footnotetext: Email: alexandre.boumezoued@milliman.com55footnotetext: Milliman, 14 Rue Pergolèse, 75016 Paris, France.

Keywords: Libor Market Model; Stochastic Volatility; Displaced Diffusion; Swaption pricing; Model calibration; Edgeworth expansions; Gram-Charlier expansions.

1 Introduction

Our work is motivated by the need in the insurance and banking industry to perform repeated calibrations of financial models. So-called market consistent forecasts are notably required for a variety of topics faced by insurance companies, such as the projection of insurance assets and liabilities, the computation of the Solvency Capital requirement through Nested Simulations, see DEVINEAU2009 and BAUER2012, the implementation of intensive recalibration process within a Least Squares Monte Carlo framework, see DEVINEAU2013, as well as for the hedging of Variable Annuities and the computation of trading grids. Among the financial models required, those dedicated to interest rates have reached a significant complexity within the insurance market practice compared to those dedicated to other financial drivers, such as stocks and inflation. Our general purpose relates to the improvement of the calibration procedure of the so-called LIBOR Market Model with Stochastic Volatility, denoted SV-LMM, which is now widely used as it has proven its ability to reproduce volatility smile and fit market prices in a satisfactory way. Additionally, in a very low interest rate regime, the use of a displacement coefficient allowing to forecast interest rates in the negative region is becoming a market standard, leading us to study the Displaced Diffusion SV-LMM, denoted DD-SV-LMM in what follows. In this context this is crucial to get fast calibration procedures, especially when the displacement coefficient itself is included in the calibration process, as such studies require to perform intensive recalibration of this coefficient in order to avoid optimization pitfalls.

Starting from the LIBOR Market Model, JOSHI2003 extended this framework to both stochastic volatility and displaced diffusion, whereas WU2006 proposed a version of the stochastic volatility component which is now widely used; on this basis they provided several analytical results such as integral-based formulas for caplets and swaptions. Several other versions of the SV-LMM have been developed in the literature, whose differences mainly lie in the way of modelling the stochastic volatility component and the scope of instruments to be addressed; for other versions of the model, we refer to references in BRIGO2007.

Due to the need for intensive repeated calibration of the model, there is a huge interest in overcoming the not-so-fast and sometimes unstable existing calibration procedures. In WU2006, pricing under the SV-LMM is performed based on both the classical HESTON1993 numerical integration method and the famous Fast Fourier Transform (FFT) approach of CARR1999, which has become a standard for option valuation for models with known characteristic function, as it is particularly the case for affine diffusion processes. Although the FFT method leads to a slight reduction (29%) in computational time compared to the Heston approach in the specific WU2006 pricing example on a strike grid (see Table 4), both methods rely on numerical integration in the complex field, which is known to embed some numerical instabilities, as already highlighted in KAHL2005 and ALBRECHER2006 on the example of the Heston model. Additionally, the numerical cost shown by both methods makes repeated calibration procedures out of reach in a reasonable operational time.

To address this issue and propose a more efficient calibration method for the DD-SV-LMM, the aim of this paper is to bring together two research areas; first, the results regarding the SV-LMM since the work of WU2006, especially on the moment generating function, and second the use of density distribution approximation based on Edgeworth and Gram-Charlier expansions. Although an analytical expression for the moment generating function does not exist for the SV-LMM in the general setting, for piecewise constant input parameters however (which are natural in the general practice), recursive closed-forms can be given, see WU2006, Proposition 4.1. This is our purpose to take advantage of this analytical tractability and implement expansions avoiding as much as possible numerical derivation and integration. This way, we perform the analytical derivation of moments up to fourth order, based on an analytical differentiation of the moment generating function. This allows us to fully exploit the potential of Edgeworth and Gram-Charlier expansions, which can be seen as an adjustment of the Bachelier model for skewness and kurtosis.

In this spirit, several contributions proposed to adjust models as primarily the Black-Scholes one for non-normal skewness and kurtosis, to overcome the well known strike price biases embedded in the standard Black-Scholes formula for away-from-the-money options. JARROW1982 derived an option pricing formula based on an Edgeworth expansion of the log-normal distribution, whereas later on, CORRADO1996 used a Gram-Charlier expansion of the normal density of log-returns in the same modelling framework. Both papers provided convincing numerical results. In our setting, we develop expansions based on the reference normal distribution; this has the advantage of providing an extension of the Bachelier model, which is our natural reference setting allowing to quote derivative instruments, as caps and swaptions, in a negative rates context; currently, short term swaption volatilities can no longer be computed in the alternative log-normal framework proposed by the Black model. Also, POTTERS1998 worked in the framework of Edgeworth expansions. They used a normal density adjusted for skewness and kurtosis, derived an analytical approximation of the volatility as a function of the cumulants, then directly fitted to the observed volatility smile (instead of prices), in an analysis dedicated to stock derivatives; this contribution is a key source of inspiration for our present study. More recent references addressed the use and/or analysis of expansions for financial models in different contexts, see e.g. SCHLOGL2013, CHATEAU2014 and HESTON2016.

By bringing together these two fields, our approach avoids the complexity and robustness issues of numerical integration, while shortening the calibration process in a significant way. A key step in the analytical tractability is the explicit derivation of moments up to fourth order which are used thereafter in the Gram-Charlier and Edgeworth expansions. Under our expansion regime, we moreover derive smile formulas relating the volatility to the moneyness. In addition to a faster calibration procedure, this therefore provides additional insights on key features on the volatility smile based on interpretable parameters. As a main conclusion, our numerical results show a 98% reduction in computational time in the DD-SV-LMM calibration process compared to the classical Fast Fourier Transform.

Our paper is structured as follows. In Section 2, we briefly sketch the swap rate dynamics underlying the DD-SV-LMM, and then proceed with the study of the moment generating function. Section 3 establishes the swaption pricing formula based on Gram-Charlier and Edgeworth expansions, and provides the related smile formulas. Finally, Section 4 details our numerical results assessing the efficiency of the proposed calibration method in comparison to the classical Heston approach. The paper ends with some concluding remarks.

2 Swap rate distribution under the DD-SV-LMM

In this section, we briefly sketch the swap rate dynamics under the Libor Market Model with Stochastic Volatility and Displaced Diffusion, denoted DD-SV-LMM in what follows. We then present the approximate swap rate dynamics under our normal volatilities framework and displaced diffusion setting, based on an adaptation of the freezing technique. Finally, we detail the set of key results on the moment generating function which will be useful to derive the analytical approximations in the next Section 3.

Although these derivations are new in this context, we omit the steps of the reasoning which are analogous to those presented in WU2006, and we refer the reader to this paper for more details.

2.1 The DD-SV-LMM framework

Let $P(t,T)$ be the zero-coupon bond maturing at time $T>t$ with par value 1. Let us introduce $F_{j}(t),\;j=1,...,M$ the value at time $t$ of the simply compounded forward rate for a period $[T_{j},T_{j+1}]$ with length $\Delta T_{j}=T_{j+1}-T_{j}$. The forward rates and zero-coupon bond prices are related through

[TABLE]

In a very low interest rate regime, the use of a displacement coefficient allows for modelling and forecasting interest rates in the negative region. Let us introduce the displacement coefficient $\delta\geq 0$, also called shift, and the $\delta$-displaced forward rate $F_{j}(t)+\delta$. The displacement coefficient $\delta$ accounts for possibly negative forward rate $F_{j}(t)$, while allowing for a log-modelling of $F_{j}(t)+\delta$. Let us introduce the forward measure ${\mathbb{Q}}^{j+1}$ associated with the numeraire $P(t,T_{j+1})$; under ${\mathbb{Q}}^{j+1}$, the displaced forward rate follows the dynamics:

[TABLE]

where the inner product ’$\cdot$’ involves a volatility vector $\zeta_{j}(t)$ and a multi-dimensional Brownian motion under ${\mathbb{Q}}^{j+1}$, denoted $Z^{j+1}$. In what follows, we denote by $m(t)=\inf\{j\geq 1:t\leq T_{j}\}$ the first forward rate that has not expired by $t$. In the model, the stochastic volatility component is specified as $\zeta_{j}(t)=\sqrt{V(t)}\gamma_{j}(t)$, where $\gamma_{j}(t)$ is a deterministic vector and $V(t)$ lies in the family of Cox-Ingersoll-Ross processes under the spot Libor measure ${\mathbb{Q}}$ associated with the numeraire $B(t)=\frac{P(t,T_{m(t)})}{\prod_{i=0}^{m(t)-1}P(T_{i},T_{i+1})}$ (sometimes assimilated to the risk neutral measure):

[TABLE]

whose Feller condition $2\kappa\theta>\epsilon^{2}$ ensures that the process has a stationary distribution and remains strictly positive. From Equation (1), it is possible to derive the stochastic dynamics of displaced forward rates under the reference risk neutral measure as, for $t\leq T_{j}$,

[TABLE]

with

[TABLE]

where $Z$ is a multi-dimensional Brownian motion under ${\mathbb{Q}}$, and correlation between $Z$ and $W$ is specified through

[TABLE]

2.2 Swap rate dynamics

Although our study can be adapted to the calibration of the model on caplets without restriction, we rather consider in this paper the calibration of the DD-SV-LMM on swaption volatilities, as it allows us to take into account correlations between forward rates. To do so, we revisit the swaption pricing as proposed in WU2006, here adapted to our setting. The swap forward rate at time $t$ for the period from $T_{m}$ to $T_{n}$ writes

[TABLE]

where $B^{S}(t)=\sum_{j=m}^{n-1}\Delta T_{j}P(t,T_{j+1})$ is the annuity of the swap (which strictly depends on $m$ and $n$ although we omit the notation for simplicity). As a numeraire, $B^{S}(t)$ defines the forward swap measure ${\mathbb{Q}}^{S}$; then the price at time zero of the payer swaption contract with strike $K$ is given by the following expectation under ${\mathbb{Q}}^{S}$:

[TABLE]

Using weights $\alpha_{j}(t)=\frac{\Delta T_{j}P(t,T_{j+1})}{B^{S}(t)}$, the swap rate can be rewritten as $R_{m,n}(t)=\sum_{j=m}^{n-1}\alpha_{j}(t)F_{j}(t)$. To value the swaption, the dynamics under ${\mathbb{Q}}^{S}$ can then be derived as follows (see WU2006, Eq. 3.3):

[TABLE]

with ${\tilde{\xi}}^{S}(t)=1+\frac{\epsilon}{\kappa}\sum_{j=m}^{n-1}\alpha_{j}(t)\sum_{k=m(t)}^{j}\frac{\Delta T_{k}(F_{k}(t)+\delta)\rho_{k}(t)\left\|\gamma_{k}(t)\right\|}{1+\Delta T_{k}F_{k}(t)}.$ The differential of the swap rate with respect to $F_{j}$ is moreover given by

[TABLE]

At this point, one faces the complexity of the dynamics, as in particular the forward rates are involved in the drift of the stochastic volatility process $V$. Analogously to ANDERSEN2000, we will proceed with the freezing technique which relies on the assumption of low variability of frozen coefficients.

Moreover, as we aim to model the swap volatility in a normal framework, we here adapt the freezing technique by fixing

[TABLE]

instead of $\frac{\partial R_{m,n}(0)}{\partial F_{j}}\frac{F_{j}(0)+\delta}{R_{m,n}(0)}$ as it would be the case in a log-normal framework. This way, we are able to approximate the swap rate dynamics as follows:

[TABLE]

where ${\tilde{\xi}}^{S}_{0}(t)=1+\frac{\epsilon}{\kappa}\sum_{j=m}^{n-1}\alpha_{j}(0)\sum_{k=m(t)}^{j}\frac{\Delta T_{k}(F_{k}(0)+\delta)\rho_{k}(t)\left\|\gamma_{k}(t)\right\|}{1+\Delta T_{k}F_{k}(0)}$.

In our setting, we develop expansions based on the reference normal distribution; this has the advantage of providing an extension of the Bachelier model, which is our natural reference setting allowing to quote derivative instruments in a negative rates context; currently, short term swaption volatilities can no longer be computed in the alternative log-normal framework associated to the Black model.

In this slightly adapted framework, it would still be possible to perform swaption pricing under the well known method developed by HESTON1993 based on numerical integration involving the characteristic function, see e.g. Equation (2.13) in WU2006. However, such approach requires the computation of an integral in the complex field, which is known to embed some possible numerical instabilities, as already highlighted in KAHL2005 and ALBRECHER2006 on the example of the Heston model. Additionally, the computational complexity involved in the numerical scheme makes repeated calibration processes out of reach in a reasonable operational time.

To address this issue and propose a more efficient calibration method for the DD-SV-LMM, we aim at providing analytical approximations of the swap rate density distribution by means of Edgeworth and Gram-Charlier expansions, leading to an adjustment of the famous Bachelier formula for skewness and kurtosis. Before detailing our expansion approach, we recall and adapt in the next subsection useful results on the moment generating function.

Remark 1**.**

When additionally one is interested into computing prices for an extended grid of strikes, the problem can be reformulated into computing a collection of summations to which the famous Fast Fourier Transform (FFT) method by CARR1999 can be applied, see e.g. Equation (5.3) in WU2006. In our study, we use as a basis for comparison of the calibration efficiency the classical method developed by HESTON1993, as indeed we will consider a limited number of strikes for out-of-the-money swaptions. As such, benchmarking with the FFT method is out of scope of the present study, and similar comparison results must hold as the orders of magnitude of the computation speed of the FFT and the Heston methods are close, see Table 4 in WU2006, and given that our 98% reduction compared to the Heston approach is significant, see Section 4 for more details. Finally, it is worth mentioning that our pricing method and smile formulas based on Gram-Charlier and Edgeworth expansions provide analytical approximations which explicitly depend on the moneyness, therefore avoiding the need for any numerical integration, and as a consequence any use of the FFT method.

2.3 The moment generating function

We present here the analytical results regarding the moment generating function in the normal volatilities framework, in which the underlying variable to characterize is the swap forward rate itself, and in our drifted diffusion setting. Let us denote by $\psi$ the moment generating function of the state variable $R_{m,n}(T_{m})$, defined by

[TABLE]

Using the fact that the conditional expectation above is a martingale, then applying ItôÕs formula and finally identifying the drift term leads to the so-called Kolmogorov backward equation

[TABLE]

with notations

[TABLE]

and terminal condition $\psi\left(x,V,T_{m};z\right)=e^{zx}$. Let us remark that this equation differs from the one exhibited in WU2006 as in the normal volatilities framework, we directly focus on the underlying process $R_{m,n}$ instead of $\ln\left(R_{m,n}+\delta\right)$. For this reason the term $-\frac{1}{2}\lambda^{2}V\frac{\partial\psi}{\partial x}$ which would appear by applying ItôÕs lemma to the process $\ln\left(R_{m,n}+\delta\right)$ vanishes in Equation (6). Adapting HESTON1993 to our context, one gets a separable form solution, with notation $\tau=T_{m}-t$,

[TABLE]

where

[TABLE]

with boundary conditions $A(0,z)=0$, $B(0,z)=0$. Note that the term $\frac{1}{2}\lambda^{2}z^{2}$ replaces the quantity $\frac{1}{2}\lambda^{2}\left(z^{2}-z\right)$ which would appear in a log-normal volatilities framework. From HESTON1993, it is possible to get an analytical closed-form expression of $A$ and $B$ under the assumption of piece-wise constant functions $\lambda$ and $\rho$ on the grid $(\tau_{j},\tau_{j+1}]$, with notation $\tau_{j}=T_{m}-T_{m-j}$, which is relevant in practice. The following recursive backward algorithm allows to compute $A$ and $B$ solution to (8): for each $j=0,...,m-1$, with convention $T_{0}=0$,

[TABLE]

where ${\tilde{A}}_{j}$ and ${\tilde{B}}_{j}$ are detailed in Appendix LABEL:appendix_solving.

3 Swaption pricing and volatility smile derived from Gram-Charlier and Edgeworth expansions

We present in this section analytical approximations for swaption prices in the DD-SV-LMM framework, allowing to extend the standard Bachelier formula to account for option smiles. The closed-forms rely on Gram-Charlier and Edgeworth expansions at fourth order, which adjust a reference Gaussian distribution by considering skewness and kurtosis. In a first step, we recall some background on Gram-Charlier and Edgeworth expansions, and discuss their main common features and differences. We then derive analytical approximations for swaption prices based on these expansions, and the closed-form derivation of moments of the swap rate up to fourth order. Finally, we develop smile formulas relating implied volatilities to the moneyness level.

3.1 Gram-Charlier and Edgeworth expansions

A Gram-Chalier series expansion (type A) of some density $f$ is defined as

[TABLE]

where $\varphi$ is the standard normal density, the $(c_{n})$ are constants related to $f$, and the $(H_{n})$ are the Hermite polynomials such that $H_{0}(z)=1$ and for $n\geq 1$,

[TABLE]

Note that for $i\neq j$, the Hermite polynomials $H_{i}$ and $H_{j}$ are orthogonal for the inner product in $L^{2}({\mathbb{R}})$ defined as ${\langle}F,G{\rangle}=\int_{\mathbb{R}}F(z)G(z)\varphi(z){\mathrm{d}}z$, allowing to identify the coefficients $(c_{n})$ which are used in what follows; the proof is left to the reader.

We consider in our study an expansion up to fourth order so as to adjust the reference density for the skewness and kurtosis of the distribution to be estimated, analogously to e.g. CORRADO1996 and NECULA2016 where Gram-Charlier series are used to adjust the Black-Scholes formula for equity option prices. Starting from a random variable $X$ of interest, with standard deviation $\nu$, we consider the density $f$ of the standardized random variable

[TABLE]

Denoting the third and fourth order moments of $Z$ by $\mu_{3}={\mathbb{E}}[Z^{3}]$ and $\mu_{4}={\mathbb{E}}[Z^{4}]$ respectively, the fourth order Gram-Charlier approximation, which we denote $g_{1}$, can be written as

[TABLE]

On the other hand, several papers rather focused on Edgeworth-type expansions, see e.g. BALIEIRO2004. Originally, and in most applications, Edgeworth expansions are used to provide an approximation of a standardized sum $S_{n}=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}X_{j}$, with $(X_{j})$ a sequence of i.i.d. standardized random variables with third and fourth order moments denoted $\gamma_{3}={\mathbb{E}}[X_{1}^{3}]$ and $\gamma_{4}={\mathbb{E}}[X_{1}^{4}]$ respectively. According to the central limit theorem, the quantity ${\mathbb{P}}(S_{n}\leq x)$ converges towards the cumulative distribution function $\Phi(x)=\int_{-\infty}^{x}\varphi(y){\mathrm{d}}y$ of the standard normal distribution. The aim of the Edgeworth expansion is to characterize the distribution of $S_{n}$ for large $n$. The second order Edgeworth expansion is often considered, which leads to the following approximation using Hermite polynomials introduced in Equation (9):

[TABLE]

Note that the term $n$ doesnÕt appear in many sudies which aim to derive pricing formulas. This issue is discussed in BALIEIRO2004, where the authors indicate that the term $n$ is incorporated to skewness and kurtosis coefficients, and leave these considerations to the reader; we propose to further detail these aspects in Appendix LABEL:appendix_edgeworth.

Let us now provide the approximation based on the single standardized random variable given in Equation (10) as

[TABLE]

Finally, after differentiation, one recovers an Edgeworth approximated density as

[TABLE]

where the density $g_{1}$ is the Gram-Charlier density introduced in (11).

3.2 Swaption pricing

The standardized random variable of interest is now

[TABLE]

with $\nu$ the standard deviation of the swap rate $R_{m,n}(T_{m})$. The price of the related swaption given in (5) now writes

[TABLE]

Let us denote indifferently $g$ the density approximation based on a Gram-Charlier or an Edgeworth expansion, as considered in (11) and (12) respectively, and still use the notations $\mu_{3}$ and $\mu_{4}$ for the third and fourth order moments of $Z$. Let us introduce the standardized moneyness $z_{k}=\frac{K-R_{m,n}(0)}{\nu}$; then the swaption price given in (14) can be approximated by

[TABLE]

Before stating our main result on swaption pricing below, let us recall that the famous Bachelier price ${\widetilde{PS}}_{0}(0,K)$ can be obtained by considering a standard normal distribution for $Z$, leading to

[TABLE]

where we recall that $\varphi$ and $\Phi$ respectively denote the standard normal density and cumulative distribution function.

Proposition 1**.**

The Gram-Charlier swaption price is given by

[TABLE]

and the Edgeworth swaption price writes

[TABLE]

According to the Newton binomial formula, the third and fourth order moments of the standardized variable $Z$ defined in Equation (13) are given as follows, for $k\in\{3,4\}$,

[TABLE]

where $\psi^{(0)}(z)\equiv\psi(R_{m,n}(0),V,0;z)$, with $V\equiv V(0)$, is the moment generating function given in analytical form in Equation (7), whose derivatives are given by

[TABLE]

where the computation of the maps $A_{m}$ and $B_{m}$ and their derivatives is detailed in Appendix LABEL:appendix_moments.

Equations (17) and (18) present the additional terms allowing to adjust the swaptions pricing Bachelier formula by taking into account the skewness and kurtosis of the swap forward rate distribution. Note that the last term of the right-hand side in (18) stems from the additional quantity which appears in the Edgeworth expansion compared to the Gram-CharlierÕs one, see Equation (12). For both formulas, one can check that any distribution with same skewness and kurtosis than the normal distribution are such that $\mu_{3}=\mu_{4}-3=0$, which makes the additional terms vanish.

Remark 2**.**

In the case of at-the-money (ATM) swaptions for which $K=R_{m,n}(0)$, the standardized moneyness $z_{K}$ is null, so that

[TABLE]

and

[TABLE]

This first shows that even for ATM swaptions, adjusted swaption prices do not match the Bachelier valuation. Furthermore, one can notice that in this case the Gram-Charlier price does not depend on the skewness of the swap rate, whereas the Edgeworth price does through the quantity $\mu_{3}^{2}$.

We now state the proof of Proposition 1.

Proof 1**.**

Let us first note that according to the property (9) on Hermite polynomials, one can get

[TABLE]

Now, let us express the Gram-Charlier density (11) into the approximated price (15), leading to

[TABLE]

where the first component reduces to the Bachelier formula (16). As for the others, it remains to compute quantities of the following form, for $j=3,4$,

[TABLE]

By integration by parts and property on Hermite polynomials, see Equation (9), the first term writes

[TABLE]

As for the second term, one gets by the property (9),

[TABLE]

This finally leads to

[TABLE]

which proves the Gram-Charlier price formula (17). The Edgeworth price formula (18) can be obtained in a similar way. The derivatives in (19) of the moment generating function can be derived by standard differentiation of (7); this is detailed in Appendix LABEL:appendix_moments.

3.3 Smile formula

In some calibration frameworks, the underlying target function to minimize in order to estimate the parameters of the DD-SV-LMM is based on volatilities instead of prices. In such a case, it may be useful to consider a smile function rather than inverting theoretical prices with a Bachelier formula. We define by smile function a closed-form expression resulting from the conversion of the Gram-Charlier or Edgeworth prices into an implied Bachelier volatility. The approach we detail hereafter to build such a smile function for swaptions instruments is an adaptation of the method proposed by BOUCHAUD2003 and BOUCHAUD2012 for stock implied volatilities.

Let us denote by $s(\nu,z_{K})$ the additive correction applied to the volatility $\nu$ in order to recover an implied Bachelier volatility, denoted $\nu(z_{K})=\nu+s(\nu,z_{K})$. Formally, based on the adjusted volatility, the Bachelier price in (16) now writes $B^{S}(0)h\left(\nu+s(\nu,z_{K})\right),$ where the function $h$ is given by

[TABLE]

The derivative of $h$ at point $\nu$ can be computed as $h^{\prime}(\nu)=\varphi(z_{K})$, which leads to the first order approximation:

[TABLE]

On the other hand, one can write the Gram-Charlier and Edgeworth prices in Equations (17) and (18), which leads to the following result.

Proposition 2**.**

The Gram-Charlier smile formula is given by

[TABLE]

and the Edgeworth smile formula writes

[TABLE]

where the third and fourth order moments $\mu_{3}$ and $\mu_{4}$ of the swap rate are mentioned in Proposition 1 and detailed in Appendix LABEL:appendix_moments.

Remark 3**.**

In the case of ATM swaptions, $K=R_{m,n}(0)$, then the standardized moneyness $z_{K}$ is null, so that the implied volatilities become

[TABLE]

This shows that the volatility of ATM swaptions $\nu_{1}(0)$ or $\nu_{2}(0)$ do not match the forward rate volatility $\nu$. Note in addition that the skewness is only involved in the Edgeworth expansion. These comments are in line with the previous Remark 2 about ATM prices.

4 Numerical results

This section details the numerical results obtained from the implementation of the expansion methods of Propositions 1 and 2, and compared to the Heston method illustrated in WU2006. We first provide an overview of the calibration setting, including the parametrization of the volatility vector, as well as the market data used. We then compare the Edgeworth and the Gram-Charlier method, and finally compare our approach to the classical Heston method based on numerical integration.

4.1 Calibration setting

In our calibration framework, we consider a piecewise constant parametrization of the volatility vector, whose value on the interval $[T_{i},T_{i+1})$ is specified as $\gamma_{j}(T_{i})=\beta_{j-i+1}g(T_{j}-T_{i})$ where $g(u)=(bu+a)e^{-cu}+d$ with non-negative constants $a,b,c\text{ and }d$, and where the $\beta_{k}$ are 2-dimensional vectors with unitary Euclidian norm. As for the correlation structure between forward rates and volatilities, we consider a constant parameter $\rho_{j}(t)=\rho$, see Equation (4). Note that the displaced coefficient $\delta$ is included in the calibration process, so that the set of parameters to be estimated is $\{a,b,c,d,\kappa,\theta,\epsilon,\rho,\delta\}$ where we recall that the parameters $\kappa,\theta\text{ and }\epsilon$ are involved in the volatility dynamics, see Equation (2).

For the purpose of illustration, the market data used for the calibration of the DD-SV-LMM are made of an average interest rate structure and swaption volatilities throughout the year 2016, for both ATM and away-from-the-money swaptions. The ATM swaptions maturities and tenors considered range into $\{1,...,10,15,20,25,30\}$. For away-from-the-money swaption volatilities, we consider the same range for maturities and focus on a 10-years reference tenor; the strikes (in bps) range into +/- $\{25,50,100,150,200\}$. In the end, this amounts to consider a set of 350 volatilities to replicate.

Finally, note that the target function to be minimized for parameter inference is computed as the sum of squared differences between market and theoretical volatilities.

4.2 Comparison between the Gram-Charlier and the Edgeworth calibration

In this part, we discuss the main differences between the Gram-Charlier and the Edgeworth expansions calibration results. As shown in Propositions 1 and 2, the Edgeworth expansion leads to an additional term compared to Gram-Charlier in the analytical approximation, this term being a function of the skewness of the swap rate distribution. Moreover, for ATM swaptions, the Edgeworth expansion still accounts for the skewness whereas it vanishes in the Gram-Charlier formulas, see Remarks 2 and 3. In practice this aspect is illustrated in Figure 1 (resp. Figure 2) by the comparison of ATM (resp. away-from-the-money) empirical volatilities for the 5-years maturity. Empirical volatilities are obtained through the following process: in a first step, we perform a calibration of the DD-SV-LMM; then, forward rates are diffused with the calibrated parameters, and we deduce empirical prices from the Monte Carlo simulations; finally, by inverting the Bachelier formula, we extract empirical volatilities. The Edgeworth approach shows a better empirical fitting accuracy of market data compared to the Gram-Charlier method. On this basis, we set the Edgeworth approach as the reference one in the next part for the comparison with the classical Heston method.

4.3 Comparison between the Edgeworth and the Heston methods

The comparison of our approach with the Heston method illustrated by WU2006 is analyzed in the light of three criteria:

(i)   In a first step, we perform a market consistency analysis assessing the fitting quality of market swaption volatilities, for both theoretical volatilities (implied by the pricing formulas) and empirical volatilities (obtained by Monte Carlo simulation). For simulation, a log-Euler scheme is taken with 5000 simulation paths, which reflects an operational standard in the insurance practice, and remains reasonable to provide satisfying convergence of Monte Carlo scenarios. Furthermore, we discuss the skew profile of theoretical volatilities (implied by the pricing formula).

(ii)   In a second step, we assess the accuracy of the Edgeworth expansion by computing theoretical volatilities based on a common set of parameters.

(iii)  In a last step, we present the gain in computational time required for calibration when using Edgeworth expansion, compared to the Heston method detailed in WU2006; as a main conclusion, our numerical results show a 98% reduction in computational time in the DD-SV-LMM calibration process compared to the classical Heston method.

For all three criteria, we present in this section the results focusing on the 5-years maturity; the results for 10-years and 20-years maturities are given in Appendix LABEL:appendix_numerical.

4.3.1 Market consistency analysis

The calibration of the DD-SV-LMM is performed on market swaption volatilities for the three following methods: the Heston method, the Edgeworth expansion applied to prices (Proposition 1) and to the volatility smile (Proposition 2), respectively called Edgeworth pricing and Edgeworth smile in the following graphs.

We report in Figure 3 the ATM empirical swaption volatilities for each method and the corresponding market swaption volatilities for the 5-years maturity. This highlights that the Heston method and both Edgeworth approaches lead to close results, providing a satisfactory Monte Carlo fitting of market data.

The theoretical (resp. empirical) volatility skew for the 5-years maturity are depicted in Figure 4 (resp. Figure 5). It can be seen that the adjustment taken by the Edgeworth expansion implies similar theoretical volatility skews as for the Heston method. Moreover theoretical and empirical results lead to close volatility skews and strongly support the ability of the Edgeworth approaches to reproduce consistently market quotes.

Note that the closeness between theoretical and empirical volatilities is more generally discussed in the next section.

4.3.2 Accuracy of approximations under a fixed set of parameters

To assess the accuracy of the approximations underlying each method, we compute for a reference set of parameters theoretical volatilities induced by the Heston method based on numerical integration, and the Edgeworth pricing and smile formula methods of Propositions 1 and 2. Then we compare these theoretical elements to Monte Carlo volatilities induced by the reference parameters. We calculate a 95% confidence interval centered on empirical volatilities and study occurrences of cases where theoretical volatilities are outside confidence intervals.

For the 5-years maturity, ATM and away-from-the-money theoretical volatilities for all methods lie in the 95% confidence interval, see Figures 6 and 7. This observation supports the robustness of the approximations in both Edgeworth pricing and smile formulas. Differences between theoretical and empirical volatilities find their origin in various reasons. On the one hand, to obtain theoretical volatilities, approximations are taken (freezing technique, and numerical integration for the Heston method, or density approximation for Edgeworth expansion). On the other hand, empirical values are biased by the sampling error, as assessed by the confidence interval, and by the log-Euler discretization scheme used for Monte Carlo simulations.

The Edgeworth pricing and smile methods lead to very similar volatility profiles both for ATM and away-from-the-money swaptions, as shown in Figure 6 and 7. Based on such analysis, the impact of the approximations involved in the Edgeworth expansion can be assessed; this appears to increase with the maturity, as depicted in Appendix LABEL:appendix_numerical for 10-years and 20-years maturities. Nevertheless, differences between theoretical and Monte Carlo swaption volatilities are small in most cases and strongly back the underlying approximations of Edgeworth expansions.

4.3.3 Gain in computational time

Numerical results of this paper have been performed under R 3.2.0, using C++ integration for key functions, in a computer with 2.6 GHz Intel Core i7 CPU. As for the comparison basis, we use a fixed budget of 2500 target function calls in the optimization routine to estimate the parameters of the DD-SV-LMM over 350 swaption volatilities, as detailed in Subsection 4.1. We report in Table LABEL:Table_CPU the CPU time in seconds needed for the calibration of the DD-SV-LMM, on a common basis of 2500 iterations budget for the optimization routine. The Edgeworth method appears much faster as it provides a 98% reduction in computational time compared to the classical Heston method.