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

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


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 [13] and [4], the implementation of intensive recalibration process within a Least Squares Monte Carlo framework, see [25], 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 it 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, [20] extended this framework to both stochastic volatility and displaced diffusion, whereas [26] 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 [6].
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 [26], pricing under the SV-LMM is performed based on both the classical [17] numerical integration method and the famous Fast Fourier Transform (FFT) approach of [7], 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 reduction in computational time compared to the Heston approach in the specific [26] 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 [21] and [1] 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 to 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 [26], 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 [26], Proposition 4.1. Our purpose is 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. [19] derived an option pricing formula based on an Edgeworth expansion of the log-normal distribution, whereas later on, [10] 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, [23] 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. [24], [8] and [18].
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; these results rely on the affine property of the approximate dynamics under study, but could be extended to the wider class of so-called polynomial processes, see [11] and [15]. 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, then proceed with the study of the moment generating function in view of moments calculation, and elaborate the link with the theory of polynomial processes. 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, the last two subsections 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, and elaborate the link with the theory of polynomial processes.
Although these derivations are new in this context, we omit the steps of the reasoning which are analogous to those presented in [26], 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 ∆T j = T j+1 − T j . The forward rates and zero-coupon bond prices are related through 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 δ ≥ 0, also called shift, and the δ-displaced forward rate F j (t) +δ. The displacement coefficient δ accounts for possibly negative forward rate F j (t), while allowing for a log-modelling of F j (t) + δ if δ > −F j (t). Let us introduce the forward measure Q j+1 associated with the numeraire P (t, T j+1 ); under Q j+1 , the displaced forward rate follows the dynamics: where the inner product '·' involves a volatility vector ζ j (t) and a multi-dimensional Brownian motion under Q j+1 , denoted by Z j+1 . In what follows, we denote by m(t) = inf{j ≥ 1 : t ≤ T j } the index of the first forward rate that has not expired by t. In the model, the stochastic volatility component is specified as where γ j (t) is a deterministic vector and V (t) lies in the family of Cox-Ingersoll-Ross processes under the spot Libor measure Q associated with the numeraire B(t) = (sometimes assimilated to the risk neutral measure): whose Feller condition 2κθ ≥ 2 ensures that the process remains strictly positive a.s.. From Equation (1), it is possible to derive the stochastic dynamics of displaced forward rates under the reference risk neutral measure as, for t ≤ T j , where Z is a multi-dimensional Brownian motion under Q. Finally, the correlation structure between forward rates is captured through the following parameters: The specification of the parameters γ j (t) and ρ j (t) considered for the calibration study will be detailed in Section 4.

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 (at least from a modelling perspective 1 ) implied correlations between forward rates. To do so, we revisit the swaption pricing as proposed in [26], here adapted to our setting. The swap forward rate at time t for the period from T m to T n writes where B S (t) = n−1 j=m ∆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 Q S ; then the price at time zero of the payer swaption contract with strike K is given by the following expectation under Q S : Using weights α j (t) = ∆Tj P (t,Tj+1) B S (t) , the swap rate can be rewritten as R m,n (t) = n−1 j=m α j (t)F j (t). To value the swaption, the dynamics under Q S can then be derived as follows (see [26], Eq. 3.3): . The differential of the swap rate with respect to F j is moreover given by 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 . Analog to [2], 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 Rm,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, by denoting by R m,n (t) the approximated swap rate (as opposed to the exact swap rate R m,n (t)): . 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 [17] based on numerical integration involving the characteristic function, see e.g. Equation (2.13) in [26]. 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 [21] and [1] 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 to 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 [7] can be applied, see e.g. Equation (5.3) in [26]. In our study, we use as a basis for comparison of the calibration efficiency the classical method developed by [17], 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 [26], 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 ψ the moment generating function of the state variable R m,n (T m ), defined by ψ (R m,n (t), V (t), t; z) = E S e zRm,n(Tm) |F t , z ∈ R.
Using the fact that the conditional expectation above is a martingale, then applying Its formula and finally identifying the drift term leads to the so-called Kolmogorov backward equation with notations Let us remark that this equation differs from the one exhibited in [26] as in the normal volatilities framework, we directly focus on the underlying process R m,n instead of ln (R m,n + δ). For this reason, the term − 1 2 λ 2 (t)V ∂ψ ∂x which would appear by applying Its lemma to the process ln (R m,n + δ) vanishes in Equation (6). Adapting [17] to our context, one gets a separable form solution, with notation where with boundary conditions A(0, z) = 0, B(0, z) = 0. Note that the term 1 2 λ 2 (T m − τ )z 2 replaces the quantity 1 2 λ 2 (T m −τ ) z 2 − z which would appear in a log-normal volatilities framework.
From [17], it is possible to get an analytical closed-form expression of A and B under the assumption of piece-wise constant functions λ and ρ on the grid (τ j , τ j+1 ], with notation τ 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, whereÃ j andB j are detailed in Appendix 5.1.
The closed-form generating function will be used in this paper to provide analytical formulas for moments computation up to order 4 based on standard differentiation. This approach therefore relies on the affine property of the approximate dynamics, which leads to tractable generating function. In the next subsection, we discuss an alternative way of representing and computing moments based on the polynomial property of the process; although this part could be omitted in a first reading, this provides simplified representations, and leaves the way for further analyses beyond affine models.

2.4.
On the polynomial property. An alternative way of computing moments for the process (R m,n (t)) t≥0 is to use the theory of polynomial processes (see [11] and [15]), which provides moments representation by means of matrix exponentials. Since all affine processes are polynomial, this theory can be elaborated in our framework.
In the original papers by [11] and [15], a theory of polynomial processes is developed for homogeneous dynamics (therefore with time-independent parameters): [11] elaborate a theory for Markov dynamics with jumps, whereas [15] focuses on diffusions, non-necessarily Markov. For such models, the expectation of any polynomial function of the process reduces to a polynomial in the initial condition with equal or lower degree, which leads to a finite system of equations for the computation of any moment. As in our context the parameters are time-dependent, we propose in the following the moment calculation as a product of matrix exponentials, which generalizes the previous work, here in an in-homogeneous setting (for piecewise constant coefficients). This provides another representation of the moments of the process, which can be used as an alternative to the direct differentiation of the generating function.

Remark 2.
Although the matrix exponential representation formulas could be used to carry out the fourth order moments expansions developed in the next section, we consider the moment calculation based on differentiation of the moment generating function to be more appropriate in our context. Indeed, as affine processes are a particular class of polynomial processes, for which the generating function is analytically known in our case, it appears that in this specific context the formulas we provide in terms of moments are fully explicit, whereas the matrix exponential (in relatively high dimension) is not 'closed-form' in comparison. From a coherence perspective, we may consider that the use of this theory provides some interest for processes which are polynomial but not affine; this leaves the way for further studies on this aspect. To this extent, this part can be omitted in a first reading.
Remark 3. From a numerical perspective, the matrix exponentiation can be tackled by applying some numerical approximations by classical algorithms (such as the 'expm' package in R); note that in our case this has lead to similar results compared to those depicted in Section 4, which are based on a more explicit derivation of the four first moments.
We consider the time-inhomogeneous Markov process X = (R m,n (t), V t ) t≥0 ; the action of the infinitesimal generator associated to the process X on (τ j , τ j+1 ] can be written as: Let us denote by P k (R 2 ) the set of polynomials on R 2 , of degree at most k; one can check that the criterion G j (P k (R 2 )) ⊂ P k (R 2 ) is satisfied for k ∈ N, which is at the basis for the characterization of polynomial processes. Let us focus on a time t in the interval (τ j , τ j+1 ]. Under the assumption of piecewise constant coefficients, the dynamics of X is homogeneous Markov between τ j and τ j+1 , with semi-group denoted by P j , and according to [11], the restriction of the semigroup P j to P k (R 2 ), k ∈ N can be written as: Finally, we detail below the moment representation in terms of matrix exponential in our context. Let us denote by B(r, v) = (b 1 (r, v), ..., b N (r, v)) the canonical basis of P k (R 2 ), and consider t ∈ (τ j , τ j+1 ]; two successive conditionings at τ j and τ j−1 and the use of the Markov property for the processX = (R m,n (t), V t , t) t≥0 leads to with, according to the generator representation (see also Theorem 3.1. in [15]), This shows that the moment of order k can be computed as which extends by recursion to the following final formula: where the product is taken in increasing order -note that the order matters as in general matrix exponentials may not commute.
k can by definition be calculated by analyzing the action of the generator G j on the canonical basis. In the case k = 4, which is our framework of interest as it involves skewness and kurtosis, the canonical basis writes k can be expressed as follows: , see those related to Equation (6), and the subscript j means that we are considering the value of timedependent coefficients on the time interval (τ j , τ j+1 ] (for instance,ρ j :=ρ(τ j+1 )).
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.

Gram-Charlier and Edgeworth expansions.
A Gram-Charlier series expansion (type A) of some density f is defined as where ϕ 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 ≥ 1, Note that for i = j, the Hermite polynomials H i and H j are orthogonal for the inner product in L 2 (R) defined as F, G = R F (z)G(z)ϕ(z)dz (see e.g. [14]), allowing to identify the coefficients (c n ) as where µ k denotes the moment of order k of density f . 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. [10] and [22] 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 ν, we consider the density f of the standardized random variable The fourth order Gram-Charlier approximation, which we denote g 1 , can be written as On the other hand, several papers rather focused on Edgeworth-type expansions, see e.g. [3]. Originally, and in most applications, Edgeworth expansions are used to provide an approximation of a standardized sum S n = 1 √ n n j=1 X j , with (X j ) being a sequence of i.i.d. standardized random variables with third and fourth order moments denoted by γ 3 = E[X 3 1 ] and γ 4 = E[X 4 1 ], respectively. According to the central limit theorem, the quantity P(S n ≤ x) converges towards the cumulative distribution function Φ(x) = x −∞ ϕ(y)dy 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 (10): Note that the term n doesnt appear in many studies which aim to derive pricing formulas. This issue is discussed in [3], 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 5.2.
Let us now provide the approximation based on the single standardized random variable given in Equation (12) as Finally, after differentiation, one recovers an Edgeworth approximated density as where the density g 1 is the Gram-Charlier density introduced in (13).

Swaption pricing. The standardized random variable of interest is now
with ν the standard deviation of the swap rate R m,n (T m ). The price of the related swaption given in (5) now writes Let us denote indifferently g the density approximation based on a Gram-Charlier or an Edgeworth expansion, as considered in (13) and (14), respectively, and still use the notations µ 3 and µ 4 for the third and fourth order moments of Z. Let us introduce the standardized moneyness z k = K−Rm,n(0) ν ; then the swaption price given in (16) can be approximated by The next Proposition details our result at first order; note that the general series representation is discussed at the end of this sub-section.
Before stating our main result on swaption pricing below, let us recall that the famous Bachelier price P S 0 (0, K) can be obtained by considering a standard normal distribution for Z, leading to where we recall that ϕ and Φ denote the standard normal density and cumulative distribution function, respectively.
Proposition 1. The Gram-Charlier swaption price is given by and the Edgeworth swaption price writes According to the Newton binomial formula, the third and fourth order moments of the standardized variable Z defined in Equation (15) are given as follows, for k ∈ {3, 4}, , is the moment generating function given in analytical form in Equation (7), whose derivatives are given by Equations (19) and (20) 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 (20) stems from the additional quantity which appears in the Edgeworth expansion compared to the Gram-Charliers one, see Equation (14). For both formulas, one can check that any distribution with the same skewness and kurtosis than the normal distribution are such that µ 3 = µ 4 − 3 = 0, which makes the additional terms vanish.

Remark 4.
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 24 .
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 µ 2 3 . We now state the proof of Proposition 1.
Proof. Let us first note that according to the property (10) on Hermite polynomials, one can get Now, let us express the Gram-Charlier density (13) into the approximated price (17), leading to where the first component reduces to the Bachelier formula (18). As for the others, it remains to compute quantities of the following form, for j = 3, 4, By integration by parts and property on Hermite polynomials, see Equation (10), the first term writes As for the second term, one gets by the property (10), This finally leads to which proves the Gram-Charlier price formula (19). The Edgeworth price formula (20) can be obtained in a similar way. The derivatives in (21) of the moment generating function can be derived by standard differentiation of (7); this is detailed in Appendix 5.3.

Remark 5.
In order to provide the reader with the untruncated Gram-Charlier series, it is worth mentioning that combining Equations (13 and 22) leads to where the c n are given in Equation (11) (we recall that c 1 = c 2 = 0 since the random variable is standardized). This formula is analogous to [24] where option pricing is performed using the Gram-Charlier framework. As for the density approximation obtained with an Edgeworth expansion, this can formally be written as where D is the differential operator and P j is a polynomial function, of degree 3j. However, to our knowledge no closed-form representation exists for the polynomial function P j (see [16]).
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 [5] and [12] for stock implied volatilities.
Let us denote by s(ν, z K ) the additive correction applied to the volatility ν in order to recover an implied Bachelier volatility, denoted by ν(z K ) = ν + s(ν, z K ). Formally, based on the adjusted volatility, the Bachelier price in (18) now writes B S (0)h (ν + s(ν, z K )) , where the function h is given by The derivative of h at point ν can be computed as h (ν) = ϕ(z K ), which leads to the first order approximation: On the other hand, one can write the Gram-Charlier and Edgeworth prices in Equations (19) and (20), which leads to the following result.
Proposition 2. The Gram-Charlier smile formula is given by and the Edgeworth smile formula writes where the third and fourth order moments µ 3 and µ 4 of the swap rate are mentioned in Proposition 1 and detailed in Appendix 5.3.

Remark 6.
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 This shows that the volatility of ATM swaptions ν 1 (0) or ν 2 (0) do not match the forward rate volatility ν. Note in addition that the skewness is only involved in the Edgeworth expansion. These comments are in line with the previous Remark 4 about ATM prices.
4. Numerical results. This section details model calibration based on market data, using the expansion methods of Propositions 1 and 2, and compares the method proposed in this paper with the Heston method illustrated in [26].

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 γ j (T i ) = β j−i+1 g(T j − T i ) where g(u) = (bu + a)e −cu + d with non-negative constants a, b, c and d, and where the β k are 2-dimensional vectors with unitary Euclidian norm. To ensure a reasonable calibration exercise, the β k parameter values have preliminarily been captured based on a Principal Component Analysis of series of log-ratios of shifted forward rates; they are then used as support vectors in the calibration process.
The correlation structure between rates and volatilities, see Equations (2) and (3), is specified as follows: the two-dimensional Brownian motion Z t = Z (1) t , Z (2) t being given, see Equation (4), an input correlation parameter ρ is considered to represent the Brownian motion W as: where Z (1) , Z (2) and Z (3) are uni-dimensional independent Brownian motions. In this setting, the correlation structure between forward rates as defined in Equation (4) can be derived as .
Note that the displaced coefficient δ is included in the calibration process, so that the set of parameters to be estimated is {a, b, c, d, κ, θ, , ρ, δ} where we recall that the parameters κ, θ and 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 336 volatilities to replicate.
Finally, note that the target function to be minimized for parameter inference is computed as the sum of squared differences (L 2 -norm) between market and theoretical volatilities -summary statistics as the sum of squared differences are displayed for each plot. Table 1, which corresponds to the sum of squared differences between Edgeworth or Gram-Charlier volatilities (which we refer to as theoretical), compared to market volatilities. As a purpose of comparison, the calibration is also performed based on the Heston method as illustrated in [26]. Several methods are considered for this comparison: the Gram-Charlier 'pricing' (Proposition 1), the Edgeworth 'pricing' (Proposition 1 again) and the Edgeworth 'smile' approach (see Proposition 2).

Summary statistics (Market vs Theoretical). The calibration target function value is given in
As expected, these target function values show that the improvement in terms of computational time (see the detailed discussion in Subsection 4.3) comes at the cost of a lower fitting accuracy to market data by both Edgeworth or Gram-Charlier methods, although the orders of magnitude of the sum of squared residuals are in
the same range, which seems reasonable from practical standards. Also, it appears that the Edgeworth approach shows a better fitting of the market volatility surface compared to the Gram-Charlier method. As shown in Propositions 1 and 2, the Edgeworth expansion formulas 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 4 and 6. These considerations are in line with the numerical results obtained. On this basis, we set the Edgeworth approach as the reference one in the next part for the comparison with the classical Heston method. Note again that these target function values represent the discrepancy between the market volatilities and those implied by the fourth order expansion formulas (Proposition 1 and Proposition 2), which we refer to as theoretical volatilities. For simulation, a log-Euler scheme is taken with 100,000 simulation paths, which lies above the operational standard in the insurance practice, and remains reasonable to provide satisfying convergence of Monte Carlo scenarios.
We report in Figure 1 the ATM Monte Carlo 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 Monte Carlo volatility skew for the 5-years maturity is depicted in Figure 2. It can be seen that the adjustment taken by the Edgeworth expansion implies similar theoretical volatility skews as for the Heston method, which supports the ability of the Edgeworth approaches to reproduce consistently market quotes. Additional results for the 10-years and 20-years maturity are given in Appendix 5.4.
Note again that the Monte Carlo simulations allow to recover 'true' model volatilities for the calibrated parameters (up to the discretization scheme and sampling error), as opposed to the theoretical volatilities obtained in Proposition 2 or through the Heston-type analysis, which involve at least the freezing techniques and/or Edgeworth-type expansions.  . In order to refine the assessment of the Edgeworth expansion accuracy, we compute for a reference set of parameters (those obtained from the Edgeworth smile calibration) the so-called 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 to Monte Carlo volatilities based on the same reference parameters -for those, a 95% confidence interval is depicted in order to gauge the reliability of the Monte Carlo volatility estimates. 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 3 and 4. 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 5.4 for 10-years and 20-years maturities. Nevertheless, differences between theoretical and Monte Carlo swaption volatilities are small in most cases.
It is worth mentioning again that differences between theoretical and Monte Carlo 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, Monte Carlo 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.   We report in Table 2 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.

Method
CPU Time (seconds) Heston 425.1 Edgeworth 8.2 Table 2. CPU time required for calibration using a 2500 optimization iterations budget This result can be mainly explained by the fact that the Heston method involves numerical integration and requires to work in the complex field, whereas the Edgeworth expansion approach takes advantage of the analytical form of swap rate moments up to fourth order, without any numerical differentiation. Furthermore, in the Edgeworth case, the derivatives of the moment generating function are only evaluated for z = 0, leading to simplified calculations.
Of course, let us emphasize again that this significant improvement in terms of numerical efficiency comes at the cost of a lower fitting accuracy compared to the Heston method, as shown in Table 1.
Finally, note that a calibration consisting in minimizing differences between market volatilities and those given by the Edgeworth smile formula is even simpler than the Edgeworth pricing method, as it doesnt require numerical inversion of the Bachelier formula during the calibration process.
The gain in speed with Edgeworth expansion enables fast recalibrations of the DD-SV-LMM and can be useful in a variety of topics faced by insurance companies, such as the computation of the Solvency Capital Requirement through Nested Simulations, see [13] and [4], the implementation of intensive recalibration process within a Least Squares Monte Carlo framework, see [25], as well as for Variable Annuities hedging and the computation of trading grids. As a matter of fact, the necessity of multiple repeated calibrations for stress-test scenarios involves a rising need for faster calibration processes. For this reason, the Edgeworth pricing and the related smile formula seem to be particularly efficient methods in an operational context. Concluding remarks. In this paper, we illustrated the efficiency of using Edgeworth and Gram-Charlier expansions applied to 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 [26], 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. The numerical results illustrated in this paper strongly back the approximations involved in the Edgeworth expansion methods, while providing satisfactory results for the fitting of market swaption volatilities. 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. It is worth mentioning again that our method works on the set of real numbers, making it much more simple and stable compared to the classical Heston approach using numerical integration in the complex field. As for further research, our method can be extended to any (even) order beyond four, so as to refine the fitting accuracy while keeping the advantage of an efficient computational approach. This could be achieved by the computation of more analytical derivatives, and a deeper understanding of the definition domain of higher order polynomials.

5.
Appendices. For simplicity of notation, the time dependency of most functionals is omitted in the following subsections. 5.1. Solving the moment generating function. Let us consider the separable form of the moment generating function introduced in Equation (7): ψ (x, V, t; z) = e A(τ,z)+B(τ,z)V +zx , with τ = T m − t. Then the first order derivatives of ψ can be computed as and second order derivatives as The Kolmogorov equation in (6) then becomes By identification, this leads to the partial differential equation in (8).

On Edgeworth expansions.
We consider in this appendix the notations introduced in Section 3. Let us recall that Edgeworth expansions are used to approximate the cumulative distribution function of a standardized sum of random variables as . (26) Note that the term n does not appear in papers which aim to derive pricing closedform expressions. This issue is discussed in [3], where the authors indicate that the term n is incorporated to skewness and kurtosis coefficients, but these considerations are left to the reader. We propose here to further detail those aspects. Denoting µ 3 = E S 3 n and µ 4 = E S 4 n , the skewness and kurtosis of S n , one recovers that µ 3 = γ3 √ n and µ 4 = γ4+3(n−1) n . Hence Equation (26) may be rewritten: This corresponds for instance to the formula used for Edgeworth Pricing adjustments in [3]. In the framework developed in Section 3, we apply this expansion to the standardized variable Z of the swap rate defined in Equation (15): Z = Rm,n(Tm)−Rm,n(0) ν . This random variable is assumed to be (finitely) divisible, that is, there exists a (possibly large) integer n and a collection of i.i.d. random increments (X j ) such that Note that this includes the set of infinitely divisible distributions, for which the previous decomposition holds for any n, as well as stable distributions which are special cases of infinitely divisible ones. As the calculation of skewness and kurtosis based on moment generating function focuses directly on the variable Z, consequently the coefficients µ 3 and µ 4 are homogeneous to those considered in our Edgeworth expansion. For this reason, we omit the term n in our Edgeworth framework.

5.3.
Moments for the swap rate distribution. Let us denote h (k) (z) = ∂ k h ∂z k (z) for any function h, with h (0) (z) = h(z), and write We recall that, denoting V = V (0) and ψ (0) (z) = ψ (R m,n (0), V, 0; z), the moment generating function writes as Let us define the following functions of z: Let us denote u = τ j+1 − τ and 3 . Functions a, d, g j and (h k ) are here implicitly time dependent; the recursive scheme is given by

FAST CALIBRATION OF THE DD-SV-LMM 1721
Note that solely the term d (0) = a (0) 2 − λ 2 2 z 2 is different to the one considered in [26]. As the state variable is R m,n (t), this leads to the additionnal R m,n (0)z term in the moment generating function. Order 1 derivative. The first derivative of the moment generating function writes where the recursive scheme for j = 0, ..., m − 1 is given by Order 2 derivative. The second derivative of the moment generating function writes where the recursive scheme for j = 0, ..., m − 1 is given by and since a (2) = 0, h Order 3 derivative. The third derivative of the moment generating function writes where the recursive scheme for j = 0, ..., m − 1 is given by Order 4 derivative. The fourth derivative of the moment generating function writes where the recursive scheme for j = 0, ..., m − 1 is given by j (z) +Ã (4) j (τ, z), ∀τ ∈ (τ j , τ j+1 ], B (4) (τ, z) = B (4) j (z) +B (4) j (τ, z), ∀τ ∈ (τ j , τ j+1 ],

5.4.
Numerical results for maturities 10-years and 20-years. We present in this appendix the numerical results for the comparison between the Edgeworth and the Heston methods; Figures 5 to 8 focus on the 10-years maturity, whereas Figures 9 to 12 are dedicated to the 20-years maturity. For each set, the first three figures relate to the calibration process, whereas the last two figures depict the swaption volatilities under a given set of parameters. The reader is referred to Section 4 for more details.   Figure 12. Swaption volatility skews with given parameters for 20-years maturity comments and feedbacks, as well as Jean-Baptiste Garnier, Abdallah Laraisse and Julien Vedani for fruitful discussions on specific topics. The authors wish to express their gratitude to the referees for their valuable comments and advice which contributed to a major improvement of the first version of the manuscript.