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

Laurent Devineau, Pierre-Edouard Arrouy, Paul Bonnefoy, Alexandre Boumezoued

► To cite this version:

Laurent Devineau, Pierre-Edouard Arrouy, Paul Bonnefoy, Alexandre Boumezoued. Fast calibration of the Libor Market Model with Stochastic Volatility and Displaced Diffusion. 2017. hal-01521491v2

HAL Id: hal-01521491

https://hal.science/hal-01521491v2

Preprint submitted on 1 Jun 2017

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Laurent Devineau<fnref n="1" />, Pierre-Edouard Arrouy<fnref n="2" />, Paul Bonnefoy<fnref n="3" />, Alexandre Boumezoued<fnref n="4" />

Milliman R&D<fnref n="5" />

June 1, 2017

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

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

<footnote n="1">¹Email: laurent.devineau@milliman.com</footnote>

<footnote n="2">²Email: pierre-edouard.arrouy@milliman.com</footnote>

<footnote n="3">³Email: paul.bonnefoy@milliman.com</footnote>

<footnote n="4">⁴Email: alexandre.boumezoued@milliman.com</footnote>

<footnote n="5">⁵Milliman, 14 Rue Pergolèse, 75016 Paris, France.</footnote>
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 Devineau and Loisel (2009) and Bauer
et al. (2012), the implementation of intensive recalibration process within a Least
Squares Monte Carlo framework, see Vedani and Devineau (2013), 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, Joshi and Rebonato (2003) extended
this framework to both stochastic volatility and displaced diffusion, whereas Wu
and Zhang (2006) 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 Brigo and
Mercurio (2007).

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 Wu and Zhang (2006), pricing under the SV-LMM is performed
based on both the classical Heston (1993) numerical integration method and the
famous Fast Fourier Transform (FFT) approach of Carr and Madan (1999), 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 Wu and Zhang (2006) 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 Kahl and Jäckel (2005) and Albrecher et al. (2006) 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 Wu and Zhang (2006), 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 Wu and Zhang (2006), 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. Jarrow and Rudd (1982) derived an option pricing formula based on an Edgeworth expansion of the log-normal distribution, whereas later on, Corrado and Su (1996) 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, Potters et al. (1998) 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. Schlögl (2013), Chateau (2014) and Heston and Rossi (2016).

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 Wu and Zhang (2006), 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

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

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 \ge 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:

$$ \text{for } t \le T_j, \quad dF_j(t) = (F_j(t) + \delta)\zeta_j(t) \cdot dZ_t^{j+1}, \qquad (1) $$

where the inner product '.' 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 \ge 1 : t \le 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):

$$ dV(t) = \kappa (\theta - V(t)) dt + \epsilon \sqrt{V(t)} dW_t, \qquad (2) $$

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 \le T_j$,

$$ dF_j(t) = (F_j(t) + \delta)\sqrt{V(t)}\gamma_j(t) \cdot \left( dZ_t - \sigma_{j+1}(t)\sqrt{V(t)}dt \right), \qquad (3) $$

with

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

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

$$ \rho_j(t)dt = \mathbb{E} \left[ \left( \frac{\gamma_j(t)}{\|\gamma_j(t)\|} \cdot dZ_t \right) dW_t \right]. \qquad (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 correlations between forward rates. To do so, we revisit the swaption pricing as proposed in Wu and Zhang (2006), here adapted to our setting. The swap forward rate at time $t$ for the period from $T_m$ to $T_n$ writes

$$ R_{m,n}(t) = \frac{P(t, T_m) - P(t, T_n)}{B^S(t)}, $$
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$:

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

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 Wu and Zhang (2006), Eq. 3.3):

$$dR_{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 dZ_t^S, \\ dV(t) = \kappa \left(\theta - \tilde{\xi}^S(t)V(t)\right) dt + \epsilon \sqrt{V(t)}dW_t^S,$$

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)\|\gamma_k(t)\|}{1+\Delta T_k F_k(t)}$. The differential of the swap rate with respect to $F_j$ is moreover given by

$$\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) (F_k(t) - R_{m,n}(t)).$$

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 Andersen and Andreasen (2000), 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

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

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:

$$dR_{m,n}(t) = \sqrt{V(t)} \sum_{j=m}^{n-1} w_j(0)\gamma_j(t) \cdot dZ_t^S, \quad 0 \le t < T_m, \\ dV(t) = \kappa \left(\theta - \tilde{\xi}_0^S(t)V(t)\right) dt + \epsilon \sqrt{V(t)}dW_t^S,$$

where $\tilde{\xi}_0^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(0)+\delta)\rho_k(t)\|\gamma_k(t)\|}{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 Heston (1993) based on numerical integration involving the characteristic function, see e.g. Equation (2.13) in Wu and Zhang (2006). 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 Kahl and Jäckel (2005) and Albrecher et al. (2006) 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 Carr and Madan (1999) can be applied, see e.g. Equation (5.3) in Wu and Zhang (2006). In our study, we use as a basis for comparison of the calibration efficiency the classical method developed by Heston (1993), 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 Wu and Zhang (2006), 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

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

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

$$ \frac{\partial \psi}{\partial t} + (\kappa\theta - \kappa\xi V) \frac{\partial \psi}{\partial V} + \frac{1}{2}\epsilon^2 V \frac{\partial^2 \psi}{\partial V^2} + \epsilon\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, \quad (6) $$

with notations

$$ \xi \equiv \tilde{\xi}_0^S(t), \quad \lambda \equiv \left\| \sum_{j=m}^{n-1} w_j(0) \gamma_j(t) \right\| \quad \text{and} \quad \rho = \frac{1}{\lambda} \sum_{j=m}^{n-1} w_j(0) \|\gamma_j(t)\| \rho_j(t), $$

and terminal condition $\psi(x, V, T_m; z) = e^{zx}$. Let us remark that this equation differs from the one exhibited in Wu and Zhang (2006) as in the normal volatilities framework, we directly focus on the underlying process $R_{m,n}$ instead of $\ln(R_{m,n} + \delta)$. 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(R_{m,n} + \delta)$ vanishes in Equation (6). Adapting Heston (1993) to our context, one gets a separable form solution, with notation $\tau = T_m - t$,

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

where

$$ \left\{ \begin{array}{l} \frac{\partial A}{\partial \tau} = \kappa \theta B, \\ \frac{\partial B}{\partial \tau} = \frac{1}{2} \epsilon^2 B^2 + (\epsilon \rho \lambda z - \kappa \xi) B + \frac{1}{2} \lambda^2 z^2, \end{array} \right. \quad (8) $$

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 (z^2 - z)$ which would appear in a log-normal volatilities framework. From Heston (1993), 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, \dots, m-1$, with convention $T_0 = 0$,

$$ \left\{ \begin{array}{l} A(\tau, z) = A(\tau_j, z) + \tilde{A}_j(\tau, z) \quad \forall \tau \in (\tau_j, \tau_{j+1}], \\ B(\tau, z) = B(\tau_j, z) + \tilde{B}_j(\tau, z) \quad \forall \tau \in (\tau_j, \tau_{j+1}], \end{array} \right. $$

where $\tilde{A}_j$ and $\tilde{B}_j$ are detailed in Appendix 5.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.

## 3.1 Gram-Charlier and Edgeworth expansions

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

$$f(z) = \varphi(z) \sum_{n=0}^{\infty} c_n H_n(z),$$

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 \ge 1$,

$$H_n(z)\varphi(z) = \varphi^{(n)}(z). \qquad (9)$$

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)dz$, 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. Corrado and Su (1996) and Necula et al. (2016) 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

$$Z = \frac{X - \mathbb{E}[X]}{\nu}. \qquad (10)$$

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

$$g_1(z) = \varphi(z) \left\{ 1 - \frac{\mu_3}{6} H_3(z) + \frac{\mu_4 - 3}{24} H_4(z) \right\}. \qquad (11)$$
On the other hand, several papers rather focused on Edgeworth-type expansions, see e.g. Balieiro Filho and Rosenfeld (2004). 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 \le x)$ converges towards the cumulative distribution function $\Phi(x) = \int_{-\infty}^{x} \varphi(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 (9):

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

Note that the term $n$ doesn't appear in many studies which aim to derive pricing formulas. This issue is discussed in Balieiro Filho and Rosenfeld (2004), 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 (10) as

$$ \mathbb{P}(Z \le z) \approx \Phi(z) - \frac{\mu_3}{6}\varphi(z)H_2(z) + \frac{\mu_4 - 3}{24}\varphi(z)H_3(z) + \frac{\mu_3^2}{72}\varphi(z)H_5(z). $$

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

$$ g_2(z) = g_1(z) + \varphi(z) \frac{\mu_3^2}{72} H_6(z), \qquad (12) $$

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

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

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

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

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

$$ \widetilde{PS}(0, K) = B^S(0) \int_{z_K}^{\infty} (R_{m,n}(0) + \nu z - K) g(z) dz = \nu B^S(0) \int_{z_K}^{\infty} (z - z_K) g(z) dz. \quad (15) $$

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

$$ \widetilde{PS}_0(0, K) = \nu B^S(0) \{\varphi(z_K) - z_K \Phi(-z_K)\}, \quad (16) $$

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

$$ \widetilde{PS}_1(0, K) = \widetilde{PS}_0(0, K) + \nu B^S(0) \varphi(z_K) \left\{ \frac{\mu_3}{6} z_K + \frac{\mu_4 - 3}{24} (z_K^2 - 1) \right\}, \quad (17) $$

and the Edgeworth swaption price writes

$$ \widetilde{PS}_2(0, K) = \widetilde{PS}_1(0, K) + \nu B^S(0) \varphi(z_K) \frac{\mu_3^2}{72} (z_K^4 - 6z_K^2 + 3). \quad (18) $$

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\}$,

$$ \mu_k = \mathbb{E}^S[Z^k] = \frac{1}{\nu^k} \sum_{j=0}^{k} \binom{k}{j} \psi^{(j)}(0) (-R_{m,n}(0))^{k-j}, $$

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

$$ \left\{ \begin{aligned} \psi^{(1)} &= (A_m^{(1)} + B_m^{(1)}V + R_{m,n}(0)) \psi^{(0)}, \\ \psi^{(2)} &= (A_m^{(2)} + B_m^{(2)}V) \psi^{(0)} + \frac{(\psi^{(1)})^2}{\psi^{(0)}}, \\ \psi^{(3)} &= (A_m^{(3)} + B_m^{(3)}V) \psi^{(0)} + (A_m^{(2)} + B_m^{(2)}V) \psi^{(1)} + \frac{2\psi^{(2)}\psi^{(1)}}{\psi^{(0)}} - \frac{(\psi^{(1)})^3}{(\psi^{(0)})^2}, \\ \psi^{(4)} &= (A_m^{(4)} + B_m^{(4)}V) \psi^{(0)} + 2(A_m^{(3)} + B_m^{(3)}V) \psi^{(1)} + (A_m^{(2)} + B_m^{(2)}V) \psi^{(2)} \\ &\quad + \frac{2\psi^{(3)}\psi^{(1)}}{\psi^{(0)}} + \frac{2(\psi^{(2)})^2}{\psi^{(0)}} - \frac{5\psi^{(2)}(\psi^{(1)})^2}{(\psi^{(0)})^2} + \frac{2(\psi^{(1)})^4}{(\psi^{(0)})^3}, \end{aligned} \right. \quad (19) $$

where the computation of the maps $A_m$ and $B_m$ and their derivatives is detailed in Appendix 5.3.
Equations (17) and (18) present the additional terms allowing to adjust the swap- tions 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

$$ \widetilde{PS}_1(0, R_{m,n}(0)) = \frac{1}{\sqrt{2\pi}} \nu B^S(0) \left\{ 1 - \frac{\mu_4 - 3}{24} \right\}, $$

and

$$ \widetilde{PS}_2(0, R_{m,n}(0)) = \frac{1}{\sqrt{2\pi}} \nu B^S(0) \left\{ 1 - \frac{\mu_4 - 3 - \mu_3^2}{24} \right\}. $$

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

$$ \left\{ \begin{array}{l} H_1(z) = -z, \\ H_2(z) = z^2 - 1, \\ H_3(z) = -z^3 + 3z, \\ H_4(z) = z^4 - 6z^2 + 3. \end{array} \right. $$

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

$$ \widetilde{PS}_1(0, K) = \nu B^S(0) \int_{z_K}^{\infty} (z - z_K) \varphi(z) dz - \nu B^S(0) \frac{\mu_3}{6} \int_{z_K}^{\infty} (z - z_K) \varphi(z) H_3(z) dz \\ + \nu B^S(0) \frac{\mu_4 - 3}{24} \int_{z_K}^{\infty} (z - z_K) \varphi(z) H_4(z) dz, $$

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

$$ \int_{z_K}^{\infty} (z - z_K) \varphi(z) H_j(z) dz = \int_{z_K}^{\infty} z \varphi(z) H_j(z) dz - z_K \int_{z_K}^{\infty} \varphi(z) H_j(z) dz, $$
By integration by parts and property on Hermite polynomials, see Equation (9), the first term writes

$$ \int_{z_K}^{\infty} z\varphi(z)H_j(z)dz = -z_K H_{j-1}(z_K)\varphi(z_K) - \int_{z_K}^{\infty} H_{j-1}(z)\varphi(z)dz. $$

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

$$ z_K \int_{z_K}^{\infty} \varphi(z) H_j(z) dz = z_K \int_{z_K}^{\infty} \varphi^{(j)}(z) dz = -z_K \varphi^{(j-1)}(z_K) = -z_K H_{j-1}(z_K) \varphi(z_K). $$

This finally leads to

$$ \int_{z_K}^{\infty} (z - z_K) \varphi(z) H_j(z) dz = - \int_{z_K}^{\infty} H_{j-1}(z) \varphi(z) dz = H_{j-2}(z_K) \varphi(z_K), $$

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

## 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 Bouchaud and Potters (2003) and De Leo et al. (2012) 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(\nu + s(\nu, z_K))$, where the function $h$ is given by

$$ h(x) = \int_{\frac{K-R_{m,n}(0)}{x}}^{\infty} (xz - (K - R_{m,n}(0)))\varphi(z)dz. $$

The derivative of $h$ at point $\nu$ can be computed as $h'(\nu) = \varphi(z_K)$, which leads to the first order approximation:

$$ h(\nu + s(\nu, z_K)) \approx h(\nu) + s(\nu, z_K)\varphi(z_K). \quad (20) $$

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

$$ \nu_1(z_K) = \nu \left\{ 1 + \frac{\mu_3}{6} z_K + \frac{\mu_4 - 3}{24} (z_k^2 - 1) \right\}, \quad (21) $$

and the Edgeworth smile formula writes

$$ \nu_2(z_K) = \nu_1(z_K) + \nu \frac{\mu_3^2}{72} (z_K^4 - 6z_K^2 + 3), \quad (22) $$

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

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

$$ \nu_1(0) = \nu \left\{ 1 - \frac{\mu_4 - 3}{24} \right\} \quad \text{and} \quad \nu_2(0) = \nu \left\{ 1 - \frac{\mu_4 - 3 - \mu_3^2}{24} \right\}. $$

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 Wu and Zhang (2006). 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$ 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$ 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 Wu and Zhang (2006) 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).
<figure><img src="image_1.png" /><figcaption>Figure 1: ATM Monte Carlo swaption volatilities for 5-years maturity</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: Monte Carlo swaption volatility skews for 5-years maturity</figcaption></figure>

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 Wu and Zhang (2006); 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 5.4.
### 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.

<figure><img src="image_3.png" /><figcaption>Figure 3: ATM Monte Carlo swaption volatilities for 5-years maturity</figcaption></figure>

### 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.
<figure><img src="image_4.png" /><figcaption>Figure 4: Theoretical swaption volatility skews for 5-years maturity</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: Monte Carlo swaption volatility skews for 5-years maturity</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>Figure 6: ATM swaption volatilities with given parameters for 5-years maturity</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Swaption volatility skew with given parameters for 5-years maturity</figcaption></figure>

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

<table><caption>Table 1: CPU time required for calibration using a 2500 optimization iterations budget</caption><thead><tr><td>Method</td><td>CPU Time (seconds)</td></tr></thead><tbody><tr><td>Heston</td><td>425.1</td></tr><tr><td>Edgeworth</td><td>8.2</td></tr></tbody></table>

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.

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 doesn't 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 Devineau and Loisel (2009) and Bauer et al. (2012), the implementation of intensive recalibration process within a Least Squares Monte Carlo framework, see Vedani and Devineau (2013), 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 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. The numeri-
cal 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 com-
pared 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 an-
alytical derivatives, and a deeper understanding of the definition domain of higher
order polynomials.
# 5 Appendices

## 5.1 Solving the moment generating function

Let us consider the separable form of the moment generating function introduced in Equation (7): $\psi(x, V, t; z) = e^{A(\tau,z)+B(\tau,z)V+zx}$, with $\tau = T_m - t$. Then the first order derivatives of $\psi$ can be computed as

$$ \frac{\partial \psi}{\partial t} = \left[ -\frac{\partial A}{\partial \tau} - V \frac{\partial B}{\partial \tau} \right] \psi, \quad \frac{\partial \psi}{\partial V} = B \psi, \quad \frac{\partial \psi}{\partial x} = z \psi, $$

and second order derivatives as

$$ \frac{\partial^2 \psi}{\partial V^2} = B^2 \psi, \quad \frac{\partial^2 \psi}{\partial x^2} = z^2 \psi, \quad \frac{\partial^2 \psi}{\partial V \partial x} = zB\psi. $$

The Kolmogorov equation in (6) then becomes

$$ \left[ -\frac{\partial A}{\partial \tau} + \kappa \theta B \right] + V \left[ -\frac{\partial B}{\partial \tau} - \kappa \xi B + \frac{1}{2} \epsilon^2 B^2 + \epsilon \rho \lambda z B + \frac{1}{2} \lambda^2 z^2 \right] = 0. $$

By identification, this leads to the partial differential equation in (8).

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

$$ \mathbb{P}(S_n \le x) \approx \Phi(x) - \frac{\gamma_3}{6\sqrt{n}}\varphi(x)H_2(x) + \frac{\gamma_4 - 3}{24n}\varphi(x)H_3(x) + \frac{\gamma_3^2}{72n}\varphi(x)H_5(x). \quad (23) $$

Note that the term $n$ does not appear in papers which aim to derive pricing closed-form expressions. This issue is discussed in Balieiro Filho and Rosenfeld (2004), 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 $\mu_3 = \mathbb{E}[S_n^3]$ and $\mu_4 = \mathbb{E}[S_n^4]$, the skewness and kurtosis of $S_n$, one recovers that $\mu_3 = \frac{\gamma_3}{\sqrt{n}}$ and $\mu_4 = \frac{\gamma_4+3(n-1)}{n}$. Hence Equation (23) may be rewritten:

$$ \mathbb{P}(S_n \le x) \approx \Phi(x) - \frac{\mu_3}{6}\varphi(x)H_2(x) + \frac{\mu_4 - 3}{24}\varphi(x)H_3(x) + \frac{\mu_3^2}{72}\varphi(x)H_5(x). $$

This corresponds for instance to the formula used for Edgeworth Pricing adjustments in Balieiro Filho and Rosenfeld (2004). In the framework developed in Section 3, we apply this expansion to the standardized variable $Z$ of the swap rate defined in Equation (13): $Z = \frac{R_{m,n}(T_m)-R_{m,n}(0)}{\nu}$. 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

$$Z = S_n = \frac{1}{\sqrt{n}} \sum_{j=1}^{n} X_j.$$

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 $\mu_3$ and $\mu_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) = \frac{\partial^k h}{\partial z^k}(z)$ for any function $h$, with $h^{(0)}(z) = h(z)$, and write

$$A_j^{(0)}(z) = A(\tau_j, z), \\ B_j^{(0)}(z) = B(\tau_j, z).$$

We recall that, denoting $V = V(0)$ and $\psi^{(0)}(z) = \psi(R_{m,n}(0), V, 0; z)$, the moment generating function writes

$$\psi^{(0)}(z) = e^{A_m^{(0)}(z) + B_m^{(0)}(z)V + zR_{m,n}(0)}.$$

Let us define the following functions of $z$:

$$a^{(0)} = \kappa\xi - \rho\epsilon\lambda z, \\ d^{(0)} = \sqrt{(a^{(0)})^2 - \lambda^2\epsilon^2 z^2}, \\ g_j = \frac{a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}}{a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)}}.$$

Let us denote $u = \tau_{j+1} - \tau$ and

$$h_1^{(0)} = a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}, \\ h_2^{(0)} = 1 - \exp(d^{(0)}u), \\ h_3^{(0)} = a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)}, \\ h_4^{(0)} = a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)} - (a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) \exp(d^{(0)}u) \\ = h_3^{(0)} + h_1^{(0)}(h_2^{(0)} - 1), \\ h_5^{(0)} = (a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) (1 - \exp(d^{(0)}u)) (a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)}), \\ = h_1^{(0)}h_2^{(0)}h_3^{(0)}.$$
Functions a, d, g_j and (h_k) are here implicitly time dependent; the recursive scheme is given by

$$
\begin{cases}
A^{(0)}(\tau, z) = A_j^{(0)}(z) + \tilde{A}_j^{(0)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \\
B^{(0)}(\tau, z) = B_j^{(0)}(z) + \tilde{B}_j^{(0)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}],
\end{cases}
$$

with

$$
\begin{align*}
\tilde{A}_j^{(0)} &= \frac{\kappa \theta}{\epsilon^2} \left[ (a^{(0)} + d^{(0)}) u - 2 \ln \left( \frac{1 - g_j \exp(d^{(0)}u)}{1 - g_j} \right) \right], \\
&= \frac{\kappa \theta}{\epsilon^2} \left[ (a^{(0)} + d^{(0)}) u - 2 \ln \left( - \frac{a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)} - (a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) \exp(d^{(0)}u)}{2d^{(0)}} \right) \right], \\
&= \frac{\kappa \theta}{\epsilon^2} \left[ (a^{(0)} + d^{(0)}) u - 2 \ln \left( - \frac{h_4^{(0)}}{2d^{(0)}} \right) \right],
\end{align*}
$$

$$
\begin{align*}
\tilde{B}_j^{(0)} &= \frac{1}{\epsilon^2} \frac{(a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) (1 - \exp(d^{(0)}u))}{1 - g_j \exp(d^{(0)}u)}, \\
&= \frac{1}{\epsilon^2} \frac{(a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) (1 - \exp(d^{(0)}u)) (a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)})}{(a^{(0)} - d^{(0)} - \epsilon^2 B_j^{(0)}) - (a^{(0)} + d^{(0)} - \epsilon^2 B_j^{(0)}) \exp(d^{(0)}u)}, \\
&= \frac{1}{\epsilon^2} \frac{h_5^{(0)}}{h_4^{(0)}}.
\end{align*}
$$

Note that solely the term $d^{(0)} = \sqrt{(a^{(0)})^2 - \lambda^2 \epsilon^2 z^2}$ is different to the one consid-
ered in Wu and Zhang (2006). 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

$$
\psi^{(1)} = (A_m^{(1)} + B_m^{(1)}V + R_{m,n}(0)) \psi^{(0)},
$$

where the recursive scheme for $j = 0, \dots, m-1$ is given by

$$
\begin{cases}
A_j^{(1)}(\tau, z) = A_j^{(1)}(z) + \tilde{A}_j^{(1)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \\
B_j^{(1)}(\tau, z) = B_j^{(1)}(z) + \tilde{B}_j^{(1)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}],
\end{cases}
$$

$$
\left\{
\begin{aligned}
\tilde{A}_j^{(1)} &= \frac{\kappa \theta}{\epsilon^2} \left[ (a^{(1)} + d^{(1)}) u - 2 \left( \frac{h_4^{(1)}}{h_4^{(0)}} - \frac{d^{(1)}}{d^{(0)}} \right) \right], \\
\tilde{B}_j^{(1)} &= \frac{1}{\epsilon^2} \left[ \frac{h_5^{(1)}}{h_4^{(0)}} - \frac{h_5^{(0)} h_4^{(1)}}{(h_4^{(0)})^2} \right],
\end{aligned}
\right.
$$
with

$$ 
\begin{aligned}
a^{(1)} &= -\rho\epsilon\lambda, \\
d^{(1)} &= \frac{a^{(0)}a^{(1)} - \lambda^2\epsilon^2z}{d^{(0)}},
\end{aligned}
$$

and

$$ 
\begin{aligned}
h_1^{(1)} &= a^{(1)} + d^{(1)} - \epsilon^2 B_j^{(1)}, \\
h_2^{(1)} &= -d^{(1)}u \exp(d^{(0)}u), \\
h_3^{(1)} &= a^{(1)} - d^{(1)} - \epsilon^2 B_j^{(1)}, \\
h_4^{(1)} &= h_3^{(1)} + h_1^{(1)}(h_2^{(0)} - 1) + h_1^{(0)}h_2^{(1)}, \\
h_5^{(1)} &= h_1^{(1)}h_2^{(0)}h_3^{(0)} + h_1^{(0)}h_2^{(1)}h_3^{(0)} + h_1^{(0)}h_2^{(0)}h_3^{(1)}.
\end{aligned}
$$

### Order 2 derivative

The second derivative of the moment generating function writes

$$ \psi^{(2)} = (A_m^{(2)} + B_m^{(2)}V) \psi^{(0)} + \frac{(\psi^{(1)})^2}{\psi^{(0)}}, $$

where the recursive scheme for $j = 0, ..., m-1$ is given by

$$ 
\left\{
\begin{aligned}
A_j^{(2)}(\tau, z) &= A_j^{(2)}(z) + \tilde{A}_j^{(2)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \\
B_j^{(2)}(\tau, z) &= B_j^{(2)}(z) + \tilde{B}_j^{(2)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}],
\end{aligned}
\right.
\\
\left\{
\begin{aligned}
\tilde{A}_j^{(2)} &= \frac{\kappa\theta}{\epsilon^2} \left[ d^{(2)}u - 2 \left( \frac{h_4^{(2)}}{h_4^{(0)}} - \frac{d^{(2)}}{d^{(0)}} - \frac{(h_4^{(1)})^2}{(h_4^{(0)})^2} + \frac{(d^{(1)})^2}{(d^{(0)})^2} \right) \right], \\
\tilde{B}_j^{(2)} &= \frac{1}{\epsilon^2} \left[ \frac{h_5^{(2)}}{h_4^{(0)}} - \frac{h_5^{(0)}h_4^{(2)}}{(h_4^{(0)})^2} - \frac{2h_5^{(1)}h_4^{(1)}}{(h_4^{(0)})^2} + \frac{2h_5^{(0)}(h_4^{(1)})^2}{(h_4^{(0)})^3} \right],
\end{aligned}
\right.
$$

with

$$ d^{(2)} = \frac{(a^{(1)})^2 - \lambda^2\epsilon^2 - (d^{(1)})^2}{d^{(0)}}, $$

and since $a^{(2)} = 0$,

$$ 
\begin{aligned}
h_1^{(2)} &= d^{(2)} - \epsilon^2 B_j^{(2)}, \\
h_2^{(2)} &= -d^{(2)}u \exp(d^{(0)}u) + d^{(1)}uh_2^{(1)}, \\
h_3^{(2)} &= -d^{(2)} - \epsilon^2 B_j^{(2)}, \\
h_4^{(2)} &= h_3^{(2)} + h_1^{(2)}(h_2^{(0)} - 1) + h_1^{(0)}h_2^{(2)} + 2h_1^{(1)}h_2^{(1)}, \\
h_5^{(2)} &= h_1^{(2)}h_2^{(0)}h_3^{(0)} + h_1^{(0)}h_2^{(2)}h_3^{(0)} + h_1^{(0)}h_2^{(0)}h_3^{(2)} + 2(h_1^{(1)}h_2^{(1)}h_3^{(0)} + h_1^{(1)}h_2^{(0)}h_3^{(1)} + h_1^{(0)}h_2^{(1)}h_3^{(1)}).
\end{aligned}
$$
## Order 3 derivative

The third derivative of the moment generating function writes

$$ \psi^{(3)} = (A_m^{(3)} + B_m^{(3)}V) \psi^{(0)} + (A_m^{(2)} + B_m^{(2)}V) \psi^{(1)} + \frac{2\psi^{(2)}\psi^{(1)}}{\psi^{(0)}} - \frac{(\psi^{(1)})^3}{(\psi^{(0)})^2}, $$

where the recursive scheme for $j = 0, ..., m-1$ is given by

$$ \begin{cases} A_j^{(3)}(\tau, z) = A_j^{(3)}(z) + \tilde{A}_j^{(3)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \\ B_j^{(3)}(\tau, z) = B_j^{(3)}(z) + \tilde{B}_j^{(3)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \end{cases} $$

$$ \left\{ \begin{aligned} \tilde{A}_j^{(3)} &= \frac{\kappa\theta}{\epsilon^2} \left[ d^{(3)}u - 2 \left( \frac{h_4^{(3)}}{h_4^{(0)}} - \frac{d^{(3)}}{d^{(0)}} - \frac{3h_4^{(2)}h_4^{(1)}}{(h_4^{(0)})^2} + \frac{3d^{(2)}d^{(1)}}{(d^{(0)})^2} + \frac{2(h_4^{(1)})^3}{(h_4^{(0)})^3} - \frac{2(d^{(1)})^3}{(d^{(0)})^3} \right) \right], \\ \tilde{B}_j^{(3)} &= \frac{1}{\epsilon^2} \left[ \frac{h_5^{(3)}}{h_4^{(0)}} - \frac{h_5^{(0)}h_4^{(3)}}{(h_4^{(0)})^2} - \frac{3h_5^{(2)}h_4^{(1)}}{(h_4^{(0)})^2} - \frac{3h_5^{(1)}h_4^{(2)}}{(h_4^{(0)})^2} + \frac{6h_5^{(0)}h_4^{(2)}h_4^{(1)}}{(h_4^{(0)})^3} + \frac{6h_5^{(1)}(h_4^{(1)})^2}{(h_4^{(0)})^3} \right. \\ & \qquad \left. - \frac{6h_5^{(0)}(h_4^{(1)})^3}{(h_4^{(0)})^4} \right], \end{aligned} \right. $$

with

$$ d^{(3)} = \frac{-3d^{(2)}d^{(1)}}{d^{(0)}}, $$

and

$$ h_1^{(3)} = d^{(3)} - \epsilon^2 B_j^{(3)}, $$

$$ h_2^{(3)} = -d^{(3)}u \exp(d^{(0)}u) + 2d^{(2)}uh_2^{(1)} + d^{(1)}uh_2^{(2)}, $$

$$ h_3^{(3)} = -d^{(3)} - \epsilon^2 B_j^{(3)}, $$

$$ h_4^{(3)} = h_3^{(3)} + h_1^{(3)}(h_2^{(0)} - 1) + h_1^{(0)}h_2^{(3)} + 3h_1^{(2)}h_2^{(1)} + 3h_1^{(1)}h_2^{(2)}, $$

$$ \begin{aligned} h_5^{(3)} = & h_1^{(3)}h_2^{(0)}h_3^{(0)} + h_1^{(0)}h_2^{(3)}h_3^{(0)} + h_1^{(0)}h_2^{(0)}h_3^{(3)} \\ & + 3(h_1^{(2)}h_2^{(1)}h_3^{(0)} + h_1^{(2)}h_2^{(0)}h_3^{(1)} + h_1^{(1)}h_2^{(2)}h_3^{(0)} + h_1^{(0)}h_2^{(2)}h_3^{(1)} + h_1^{(1)}h_2^{(0)}h_3^{(2)} + h_1^{(0)}h_2^{(1)}h_3^{(2)}) \\ & + 6h_1^{(1)}h_2^{(1)}h_3^{(1)}. \end{aligned} $$

## Order 4 derivative

The fourth derivative of the moment generating function writes

$$ \psi^{(4)} = (A_m^{(4)} + B_m^{(4)}V) \psi^{(0)} + 2(A_m^{(3)} + B_m^{(3)}V) \psi^{(1)} + (A_m^{(2)} + B_m^{(2)}V) \psi^{(2)} \\ + \frac{2\psi^{(3)}\psi^{(1)}}{\psi^{(0)}} + \frac{2(\psi^{(2)})^2}{\psi^{(0)}} - \frac{5\psi^{(2)}(\psi^{(1)})^2}{(\psi^{(0)})^2} + \frac{2(\psi^{(1)})^4}{(\psi^{(0)})^3}, $$
where the recursive scheme for $j = 0, ..., m-1$ is given by

$$ egin{cases} A_j^{(4)}(\tau, z) = A_j^{(4)}(z) + \tilde{A}_j^{(4)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \\ B_j^{(4)}(\tau, z) = B_j^{(4)}(z) + \tilde{B}_j^{(4)}(\tau, z), & \forall \tau \in (\tau_j, \tau_{j+1}], \end{cases} $$

$$ egin{cases} 
\tilde{A}_j^{(4)} = \frac{\kappa\theta}{\epsilon^2} \left[ d^{(4)}u - 2 \left( \frac{h_4^{(4)}}{h_4^{(0)}} - \frac{d^{(4)}}{d^{(0)}} - \frac{4h_4^{(3)}h_4^{(1)}}{(h_4^{(0)})^2} + \frac{4d^{(3)}d^{(1)}}{(d^{(0)})^2} - \frac{3(h_4^{(2)})^2}{(h_4^{(0)})^2} + \frac{3(d^{(2)})^2}{(d^{(0)})^2} \right. \\ 
\left. + \frac{12h_4^{(2)}(h_4^{(1)})^2}{(h_4^{(0)})^3} - \frac{12d^{(2)}(d^{(1)})^2}{(d^{(0)})^3} - \frac{6(h_4^{(1)})^4}{(h_4^{(0)})^4} + \frac{6(d^{(1)})^4}{(d^{(0)})^4} \right] \\ 
\tilde{B}_j^{(4)} = \frac{1}{\epsilon^2} \left[ \frac{h_5^{(4)}}{h_4^{(0)}} - \frac{h_5^{(0)}h_4^{(4)}}{(h_4^{(0)})^2} - \frac{4h_5^{(3)}h_4^{(1)}}{(h_4^{(0)})^2} - \frac{4h_5^{(1)}h_4^{(3)}}{(h_4^{(0)})^2} - \frac{6h_5^{(2)}h_4^{(2)}}{(h_4^{(0)})^2} + \frac{8h_5^{(0)}h_4^{(3)}h_4^{(1)}}{(h_4^{(0)})^3} + \frac{12h_5^{(2)}(h_4^{(1)})^2}{(h_4^{(0)})^3} \right. \\ 
\left. + \frac{24h_5^{(1)}h_4^{(2)}h_4^{(1)}}{(h_4^{(0)})^3} + \frac{6h_5^{(0)}(h_4^{(2)})^2}{(h_4^{(0)})^3} - \frac{36h_5^{(0)}h_4^{(2)}(h_4^{(1)})^2}{(h_4^{(0)})^4} - \frac{24h_5^{(1)}(h_4^{(1)})^3}{(h_4^{(0)})^4} + \frac{24h_5^{(0)}(h_4^{(1)})^4}{(h_4^{(0)})^5} \right] 
\end{cases} $$

with

$$ d^{(4)} = -\frac{3d^{(3)}d^{(1)}}{d^{(0)}} - \frac{3(d^{(2)})^2}{d^{(0)}} + \frac{3d^{(2)}(d^{(1)})^2}{(d^{(0)})^2}, $$

and

$$ h_1^{(4)} = d^{(4)} - \epsilon^2 B_j^{(4)}, $$

$$ h_2^{(4)} = -d^{(4)}u \exp(d^{(0)}u) + d^{(1)}uh_2^{(3)} + 3d^{(3)}uh_2^{(1)} + 3d^{(2)}uh_2^{(2)}, $$

$$ h_3^{(4)} = -d^{(4)} - \epsilon^2 B_j^{(4)}, $$

$$ h_4^{(4)} = h_3^{(4)} + h_1^{(4)}(h_2^{(0)} - 1) + h_1^{(0)}h_2^{(4)} + 4h_1^{(3)}h_2^{(1)} + 4h_1^{(1)}h_2^{(3)} + 6h_1^{(2)}h_2^{(2)}, $$

$$ 
\begin{aligned} 
h_5^{(4)} = & h_1^{(4)}h_2^{(0)}h_3^{(0)} + h_1^{(0)}h_2^{(4)}h_3^{(0)} + h_1^{(0)}h_2^{(0)}h_3^{(4)} \\ 
& + 4(h_1^{(3)}h_2^{(1)}h_3^{(0)} + h_1^{(3)}h_2^{(0)}h_3^{(1)} + h_1^{(1)}h_2^{(3)}h_3^{(0)} + h_1^{(0)}h_2^{(3)}h_3^{(1)} + h_1^{(1)}h_2^{(0)}h_3^{(3)} + h_1^{(0)}h_2^{(1)}h_3^{(3)}) \\ 
& + 6(h_1^{(2)}h_2^{(2)}h_3^{(0)} + h_1^{(2)}h_2^{(0)}h_3^{(2)} + h_1^{(0)}h_2^{(2)}h_3^{(2)}) + 12(h_1^{(2)}h_2^{(1)}h_3^{(1)} + h_1^{(1)}h_2^{(2)}h_3^{(1)} + h_1^{(1)}h_2^{(1)}h_3^{(2)}). 
\end{aligned} 
$$

## 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 8 to 12 focus on the 10-years maturity, whereas Figures 13 to 17 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><img src="image_8.png" /><figcaption>Figure 8: ATM Monte Carlo swaption volatilities for 10-years maturity</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9: Theoretical swaption volatility skews for 10-years maturity</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 10: Monte Carlo swaption volatility skews for 10-years maturity</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 11: ATM swaption volatilities with given parameters for 10-years maturity</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Swaption volatility skews with given parameters for 10-years maturity</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 13: ATM Monte Carlo swaption volatilities for 20-years maturity</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 14: Theoretical swaption volatility skews for 20-years maturity</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Monte Carlo swaption volatility skews for 20-years maturity</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 16: ATM swaption volatilities with given parameters for 20-years maturity</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 17: Swaption volatility skews with given parameters for 20-years maturity</figcaption></figure>
References

Albrecher, Hansjörg, Philipp Mayer, Wim Schoutens, Jurgen Tistaert. 2006. The little Heston trap. 3, 7

Andersen, Leif, Jesper Andreasen. 2000. Volatility skews and extensions of the Libor market model. *Applied Mathematical Finance* **7**(1) 1–32. 6

Balieiro Filho, Ruy Gabriel, Rogerio Rosenfeld. 2004. Testing option pricing with the Edgeworth expansion. *Physica A: Statistical Mechanics and its Applications* **344**(3) 484–490. 10, 22

Bauer, Daniel, Andreas Reuss, Daniela Singer. 2012. On the calculation of the Solvency Capital Requirement based on nested simulations. *Astin Bulletin* **42**(02) 453–499. 2, 20

Bouchaud, Jean-Philippe, Marc Potters. 2003. *Theory of financial risk and derivative pricing: from statistical physics to risk management*. Cambridge university press. 13

Brigo, Damiano, Fabio Mercurio. 2007. *Interest rate models-theory and practice: with smile, inflation and credit*. Springer Science & Business Media. 2

Carr, Peter, Dilip Madan. 1999. Option valuation using the fast Fourier transform. *Journal of computational finance* **2**(4) 61–73. 2, 7

Chateau, John-Peter D. 2014. Valuing European put options under skewness and increasing (excess) kurtosis. *Journal of Mathematical Finance* **4**(03) 160. 3

Corrado, Charles J, Tie Su. 1996. Skewness and kurtosis in S&P 500 index returns implied by option prices. *Journal of Financial Research* **19**(2) 175–192. 3, 9

De Leo, L, V Vargas, S Ciliberti, J-P Bouchaud. 2012. We’ve walked a million miles for one of these smiles. *arXiv preprint arXiv:1203.5703*. 13

Devineau, Laurent, Stéphane Loisel. 2009. Construction d'un algorithme d'accélération de la méthode des «simulations dans les simulations» pour le calcul du capital économique Solvabilité II. *Bulletin Français d'Actuariat* **10**(17) 188–221. 2, 20

Heston, Steven, Alberto G Rossi. 2016. A spanning series approach to options. *Review of Asset Pricing Studies*. 3

Heston, Steven L. 1993. A closed-form solution for options with stochastic volatility with applications to bond and currency options. *Review of financial studies* **6**(2) 327–343. 1, 2, 7, 8

Jarrow, Robert, Andrew Rudd. 1982. Approximate option valuation for arbitrary stochastic processes. *Journal of financial Economics* **10**(3) 347–369. 3

Joshi, Mark, Riccardo Rebonato. 2003. A stochastic-volatility, displaced-diffusion extension of the LIBOR market model. *Quantitative Finance* **3**(6) 458–469. 2

Kahl, Christian, Peter Jäckel. 2005. Not-so-complex logarithms in the Heston model. *Wilmott magazine* **19**(9) 94–103. 3, 7

Necula, Ciprian, Gabriel G Drimus, Walter Farkas. 2016. A general closed form option pricing formula. *Swiss Finance Institute Research Paper No. 15-53*. 9
Potters, Marc, Rama Cont, J-P Bouchaud. 1998. Financial markets as adaptive systems. *EPL (Europhysics Letters)* **41**(3) 239. 3

Schlögl, Erik. 2013. Option pricing where the underlying assets follow a Gram/Charlier density of arbitrary order. *Journal of Economic Dynamics and Control* **37**(3) 611–632. 3

Vedani, Julien, Laurent Devineau. 2013. Solvency assessment within the ORSA framework: issues and quantitative methodologies. *Bulletin Français d'Actuariat* **13**(25) 35–71. 2, 20

Wu, Lixin, Fan Zhang. 2006. LIBOR market model with stochastic volatility. *Journal of industrial and management optimization* **2**(2) 199. 1, 2, 3, 4, 5, 6, 7, 8, 14, 15, 16, 20, 24

## Acknowledgements

The authors are grateful to their colleagues at Milliman for fruitful discussions and enlightening comments, in particular Jean-Baptiste Garnier, Abdallah Laraisse, Damien Louvet and Julien Vedani.