Calibration of a Hybrid Local-Stochastic Volatility Stochastic Rates Model with a Control Variate Particle Method

TL;DR

This paper introduces a new calibration method for complex four-factor foreign-exchange models that combines particle methods with variance reduction, enabling efficient calibration on real market data.

Contribution

It presents a generic calibration technique for hybrid local-stochastic volatility models with stochastic rates, improving convergence speed using control variates and particle methods.

Findings

Effective calibration on EUR-USD data

Comparable run-time to PDE calibration of simpler models

Applicable to a broad class of hybrid LSV models

Abstract

We propose a novel and generic calibration technique for four-factor foreign-exchange hybrid local-stochastic volatility models with stochastic short rates. We build upon the particle method introduced by Guyon and Labord\`ere [Nonlinear Option Pricing, Chapter 11, Chapman and Hall, 2013] and combine it with new variance reduction techniques in order to accelerate convergence. We use control variates derived from a calibrated pure local volatility model, a two-factor Heston-type LSV model (both with deterministic rates), and the stochastic (CIR) short rates. The method can be applied to a large class of hybrid LSV models and is not restricted to our particular choice of the diffusion. The calibration procedure is performed on real-world market data for the EUR-USD currency pair and has a comparable run-time to the PDE calibration of a two-factor LSV model alone.

Figures

Tables

Table 1. Table 1 : Average error in absolute volatility (% unit) for a high-volatility stress scenario on the rate processes. N 𝑁 N is the number of particles.

$N$	$\mathrm{4\hspace{0.17em}000}$	$\mathrm{20\hspace{0.17em}000}$	$\mathrm{100\hspace{0.17em}000}$	$\mathrm{500\hspace{0.17em}000}$	$\mathrm{2\hspace{0.17em}500\hspace{0.17em}000}$
Plain particle method	$0.0855\%$	$0.0440\%$	$0.0320\%$	$0.0255\%$	$0.0160\%$
Control variate particle method	$0.0312\%$	$0.0165\%$	$0.0133\%$	$0.0133\%$	$0.0133\%$

Table 2. Table 2 : Monte Carlo NPVs, in % of N EUR subscript 𝑁 EUR N_{\text{EUR}} and N USD subscript 𝑁 USD N_{\text{USD}} respectively, for the 5Y foreign no-touch and target accrual redemption note

Model	NPV - No-touch	NPV - TARN
LV	$72.42\%$	$29.16\%$
LSV	$74.59\%$	$14.69\%$
LSV–2CIR ++	$73.60\%$	$17.50\%$

Table 3. Table 3 : The calibrated CIR parameters

CCY	$g_0$	$\kappa$	$\theta$	$\xi$
USD	$0.0001$	$0.0837$	$0.5469$	$0.0274$
EUR	$0.0001$	$0.0110$	$1.1656$	$0.0370$

Equations

$$D_t^d={\text{e}}^{-{\int }_0^t{r}_u^{d}du},D_t^f={\text{e}}^{-{\int }_0^t{r}_u^{f}du}.$$

$$\begin{array}[]{c}\begin{cases}\cfrac{dS_{t}}{S_{t}}=\left(r_{t}^{d}-r_{t}^{f}\right)\,dt+\alpha\left(S_{t},t\right)\sqrt{V_{t}}\,dW_{t}\\ r_{t}^{d}=g_{t}^{d}+h^{d}\left(t\right)\\ r_{t}^{f}=g_{t}^{f}+h^{f}\left(t\right)\\ dg_{t}^{d}=\kappa_{d}\left(\theta_{d}-g_{t}^{d}\right)\,dt+\xi_{d}\sqrt{g_{t}^{d}}\,dW_{t}^{d}\\ dg_{t}^{f}=\left(\kappa_{f}\left(\theta_{f}-g_{t}^{f}\right)-\rho_{Sf}\xi_{f}\sqrt{g_{t}^{f}}\alpha\left(S_{t},t\right)\sqrt{V_{t}}\right)\,dt+\xi_{f}\sqrt{g_{t}^{f}}\,dW_{t}^{f}\\ dV_{t}=\kappa\left(\theta-V_{t}\right)\,dt+\xi\sqrt{V_{t}}\,dW_{t}^{V},\end{cases}\end{array}$$

$$d⟨W_t,W_t^V⟩=\rho dt,d⟨W_t,W_t^d⟩={\rho }_{Sd}dt,d⟨W_t,W_t^f⟩={\rho }_{Sf}dt,d⟨W_t^d,W_t^f⟩={\rho }_{df}dt,$$

$$C\left(K,T\right)={\mathbb{E}}^{{\mathbb{Q}}^d}\left[D_T^d{\left(S_T-K\right)}^{+}\right].$$

$$\begin{array}[]{c}\begin{cases}\cfrac{dS_{t}^{2D}}{S_{t}^{2D}}=\left(\bar{r}^{d}\left(t\right)-\bar{r}^{f}\left(t\right)\right)\,dt+\alpha^{2D}\left(S_{t}^{2D},t\right)\sqrt{V_{t}^{2D}}\,dW_{t},\quad S_{0}^{2D}=S_{0},\\ dV_{t}^{2D}=\kappa\left(\theta-V_{t}^{2D}\right)\,dt+\xi\sqrt{V_{t}^{2D}}\,dW_{t}^{V},\quad V_{0}^{2D}=V_{0},\end{cases}\end{array}$$

$$\frac{d{S}_t^{LV}}{S_t^{LV}}=\left({\bar{r}}^d\left(t\right)-{\bar{r}}^f\left(t\right)\right)dt+{\sigma }_{LV}\left(S_t^{LV},t\right)d{W}_t,S_0^{LV}=S_0,$$

$$\frac{\partial C_{LV}}{\partial T}-\frac{1}{2}{\sigma }_{LV}{\left(K,T\right)}^2{K}^2\frac{{\partial }^2{C}_{LV}}{\partial K^2}+K\left({\bar{r}}^d\left(T\right)-{\bar{r}}^f\left(T\right)\right)\frac{\partial C_{LV}}{\partial K}+{\bar{r}}^f\left(T\right)C_{LV}=0.$$

$${\alpha }^2\left(K,T\right)$$

$$Q_T$$

$$\{\begin{array}{ll}{T}^{\ast }=\frac{2}{\sqrt{{\phi }^2{\zeta }^2-{\kappa }^2}}\left[\frac{\pi }{2}+\arctan \left(\frac{\kappa }{\sqrt{{\phi }^2{\zeta }^2-{\kappa }^2}}\right)\right],& \text{if }\kappa <\phi \zeta ,\\ T^{\ast }=\infty ,& \text{if }\kappa \ge \phi \zeta ,\end{array}$$

$$\frac{\partial C\left(K,T\right)}{\partial T}-\frac{1}{2}{\alpha }^2\left(K,T\right)K^2\frac{{\mathbb{E}}^{{\mathbb{Q}}^d}\left[D_T^d{V}_T|S_T=K\right]}{{\mathbb{E}}^{{\mathbb{Q}}^d}\left[D_T^d|S_T=K\right]}\frac{{\partial }^{2}C\left(K,T\right)}{\partial K^2}$$

$$+{\mathbb{E}}^{{\mathbb{Q}}^d}\left[D_T^d{r}_T^f{\left(S_T-K\right)}^{+}\right]-{\mathbb{E}}^{{\mathbb{Q}}^d}\left[D_T^d{1}_{S_T\ge K}K\left(r_T^d-r_T^f\right)\right]=0$$

$${\alpha }^{2D}\left(K,T\right)=\frac{{\sigma }_{LV}\left(K,T\right)}{\sqrt{{\mathbb{E}}^{{\mathbb{Q}}^d}\left[V_T^{2D}|S_T^{2D}=K\right]}}.$$

$${\int }_0^t{1}_{S_u\ge K}{D}_u^d{S}_u\alpha \left(S_u,u\right)\sqrt{V_u}d{W}_u$$

$$\{\begin{array}{ll}\frac{d{S}_t^i}{S_t^i}=\left(r_t^{d,i}-r_t^{f,i}\right)dt+{\hat{\alpha }}_N\left(S_t^i,t,{\left(X_t^{j,N}\right)}_{j\le N}\right)\sqrt{V_t^i}d{W}_t^i& \\ r_t^{d,i}=g_t^{d,i}+h^d\left(t\right)& \\ r_t^{f,i}=g_t^{f,i}+h^f\left(t\right)& \\ d{g}_t^{d,i}={\kappa }_d\left({\theta }_d-g_t^{d,i}\right)dt+{\xi }_d\sqrt{g_t^{d,i}}d{W}_t^{d,i}& \\ d{g}_t^{f,i}=\left({\kappa }_f\left({\theta }_f-g_t^{f,i}\right)-{\rho }_{Sf}{\xi }_f\sqrt{g_t^{f,i}}{\hat{\alpha }}_N\left(S_t^i,t,{\left(X_t^{j,N}\right)}_{j\le N}\right)\sqrt{V_t^i}\right)dt+{\xi }_f\sqrt{g_t^{f,i}}d{W}_t^{f,i}& \\ d{V}_t^i=\kappa \left(\theta -V_t^i\right)dt+\xi \sqrt{V_t^i}d{W}_t^{V,i}& \\ d{D}_t^{d,i}=-r_t^{d,i}{D}_t^{d,i}dt,& \end{array}$$

$${\hat{\alpha }}_N\left(K,T,{\left(X_t^{i,N}\right)}_{i\le N}\right)=\sqrt{\frac{{\sigma }_{LV}{\left(K,T\right)}^2}{{\hat{p}}_N\left(K,T\right)}+\frac{\hat{Q}\left(K,T\right)}{\frac{1}{2}{\hat{p}}_N\left(K,T\right)K^2\frac{{\partial }^2{C}_{LV}}{\partial K^2}}},$$

$$\hat{Q}=\frac{1}{N}\sum _{i=1}^N{Q}_T^i\text{and}{\hat{p}}_N\left(K,T\right)=\frac{{\sum }_{i=1}^N{D}_T^{d,i}{V}_T^i{\delta }_N\left(S_T^i-K\right)}{{\sum }_{i=1}^N{D}_T^{d,i}{\delta }_N\left(S_T^i-K\right)},$$

$${\hat{p}}_N^{2D}\left(K,T\right)=\sum _{i=1}^N\frac{V_T^{2D,i}{\delta }_N\left(S_T^{2D,i}-K\right)}{{\sum }_{i=1}^N{\delta }_N\left(S_T^{2D,i}-K\right)},$$

$$p^{2D}(K,T)={\mathbb{E}}^{{\mathbb{Q}}^d}\left[V_T^{2D}|S_T^{2D}=K\right],$$

$$p_N^{\ast }(K,T)={\hat{p}}_N\left(K,T\right)+\lambda \left({\hat{p}}_N^{2D}\left(K,T\right)-p^{2D}\left(K,T\right)\right).$$

$$p_N^{\ast }=\frac{1}{N}\sum _{i=1}^N{m}_i+\lambda \left(\frac{1}{N}\sum _{i=1}^N{m}_i^{2D}-p^{2D}\right)$$

$$\begin{array}[]{cc}m_{i}=\frac{D_{T}^{d,i}V_{T}^{i}\delta_{N}\left(S_{T}^{i}-K\right)}{\frac{1}{N}\sum_{i=1}^{N}D_{T}^{d,i}\delta_{N}\left(S_{T}^{i}-K\right)},&m_{i}^{2D}=\frac{V_{T}^{2D,i}\delta_{N}\left(S_{T}^{2D,i}-K\right)}{\frac{1}{N}\sum_{i=1}^{N}\delta_{N}\left(S_{T}^{2D,i}-K\right)}\end{array}\,,$$

$$\lambda =-\frac{\text{Cov}\left(m,m^{2D}\right)}{\text{Var}\left(m^{2D}\right)},$$

$${\hat{\lambda }}_N=-\frac{{\sum }_{i=1}^N\left(m_i-{\hat{p}}_N\right)\left(m_i^{2D}-{\hat{p}}_N^{2D}\right)}{{\sum }_{i=1}^N{\left(m_i^{2D}-{\hat{p}}_N^{2D}\right)}^2}.$$

$$\frac{1}{1-\text{Corr}{\left(m,m^{2D}\right)}^2}.$$

$$Q_T=\underset{=:X_{1,T}}{\underbrace{D_T^d\left(r_T^f-{\bar{r}}^f\left(T\right)\right){\left(S_T-K\right)}^{+}}}-K\underset{=:X_{2,T}}{\underbrace{D_T^d{1}_{S_T\ge K}\left[\left(r_T^d-{\bar{r}}^d\left(T\right)\right)-\left(r_T^f-{\bar{r}}^f\left(T\right)\right)\right]}},$$

$$Y_{1,T}=D_T^d{\left(S_T-K\right)}^{+},Z_{1,T}={\hat{r}}_T^f-{\bar{r}}^f\left(T\right)$$

$$Y_{2,T}=D_T^d{1}_{S_T\ge K},Z_{2,T}=\left(r_T^d-{\bar{r}}^d\left(T\right)\right)-\left({\hat{r}}_T^f-{\bar{r}}^f\left(T\right)\right)$$

$${\mathbb{E}}^{{\mathbb{Q}}^d}\left[Y_{1,T}\right]$$

$${\mathbb{E}}^{{\mathbb{Q}}^d}\left[Y_{2,T}\right]$$

$${\zeta }_1\equiv {\mathbb{E}}^{{\mathbb{Q}}^d}\left[Z_{1,T}\right]$$

Full text

1 Introduction

Efficient pricing and hedging of exotic derivatives requires a model which is rich enough to re-price accurately a range of liquidly traded market products. Calibration to vanilla options has been widely documented in the literature since the work of Dupire [20] in the context of local volatility (LV). Nowadays, the exact re-pricing of call options is a must-have standard, and Local-Stochastic Volatility (LSV) models are the state-of-the-art in many financial institutions. As discussed in Ren et al. [39], Tian et al. [45], Van der Stoep et al. [18] and Guyon and Henry-Labordère [26], LSV models improve the pricing and risk-management performance when compared to pure local volatility or pure stochastic volatility models. The local volatility component allows a perfect calibration to the market prices of vanilla options. At the same time, the stochastic volatility component already provides built-in smiles and skews which give a rough fit, so that a local volatility component – the so-called leverage function – relatively close to one suffices for a perfect calibration. Moreover, they exhibit superior dynamic properties over pure local volatility models.

We focus on a Heston-type LSV model because of the desirable properties of the Cox–Ingersoll–Ross (CIR) process for the variance, such as mean-reversion and non-negativity, and since semi-analytic formulae are available for calls and puts under Heston’s model (see [28]) and can help calibrate the Heston parameters easily. Various sophisticated calibration techniques for the local volatility component are in use in the financial industry, e.g., based on the Monte Carlo particle method in [26] or the PDE-based approach in [39].

In order to improve the pricing and hedging of foreign exchange (FX) options, we furthermore introduce stochastic domestic and foreign short interest rates into the model. Empirical results (see e.g. [48]) have confirmed that for long-dated FX products the effect of interest rate volatility can be as relevant as that of the FX rate volatility. Extensive research has been carried out in the area of option pricing with stochastic volatility and interest rates in the past few years. Van Haastrecht et al. [48] extended the model of Schöbel and Zhu [41] to currency derivatives by including stochastic interest rates, a model that benefits from analytical tractability even in a full correlation setting due to the processes being Gaussian. On the other hand, Ahlip and Rutkowski [3], Grzelak and Oosterlee [25] and Van Haastrecht and Pelsser [49] examined Heston–CIR/Vasicek hybrid models and concluded that they give rise to non-affine models even under a partial correlation structure of the driving Brownian motions and are not analytically tractable.

The resulting 4-factor model complicates the calibration routine due to the higher dimensionality, especially when PDEs are used to find the joint distribution of all factors. A few papers discuss this problem in simpler settings. Deelstra [17] and Clark [11] mainly consider 3-factor hybrid local volatility models and focus on the theoretical rather than the practical aspects of the calibration, whereas Van der Stoep et al. [47] consider an application to a 2-factor hybrid local volatility. In [26], Guyon and Henry-Labordère discuss an application of Monte Carlo-based calibration methods to a 3-factor LSV equity model with stochastic domestic rate and discrete dividends.

The model of Cox et al. [14] is popular when modeling short rates because the (square-root) CIR process admits a unique strong solution, is mean-reverting and analytically tractable. As of late, the non-negativity of the CIR process is considered to be less desirable when modeling short rates. On one hand, central banks have significantly reduced the interest rates since the 2008 financial crisis and it is now commonly accepted that interest rates need not be positive. On the other hand, if interest rates dropped too far below zero, then large amounts of money would be withdrawn from banks and government bonds, putting a severe squeeze on deposits. Hence, we model the domestic and foreign short rates using the shifted CIR (CIR

++

) process of Brigo and Mercurio [8]. The CIR

++

model allows the short rates to become negative and can fit any observed term structure exactly while preserving the analytical tractability of the original model for bonds, caps, swaptions and other basic interest rate products.

We note that the CIR process is sometimes considered difficult to simulate in practice. Moreover, as factor in the Heston model, it leads to singular probability densities for parameter settings where the variance process can hit zero (i.e., if the so-called Feller condition is violated), which cannot be handled easily in the forward Kolmogorov equation by standard numerical methods. In this paper, we address both these issues by tailored schemes, but note that the variance and interest rate processes can be exchanged without significant changes to the main framework and its benefits, e.g. by exponential Ornstein-Uhlenbeck processes for the volatility and Hull-White processes for the rates, both of which are also popular in the industry.

Based on the above considerations, we study the 4-factor hybrid LSV model defined in (2.1) below, which is a Heston-type LSV model with two shifted CIR short-rate processes. We give a rigorous proof of the calibration condition for the leverage function given in [27, Proposition 12.8] for our model specification; see also the condition given in [17] for a 4-factor LSV–2CIR++ model.

We propose a calibration approach which builds on the particle method of [26], and combines it with a novel and efficient variance reduction technique. The main control variate is the two-factor LSV model obtained by assuming that the domestic and foreign rates are deterministic in the original model. In this case, the leverage function is computed by using a deterministic PDE solver. This allows us to take advantage of the efficiency and accuracy of PDE calibration for a low-dimensional model while keeping the complexity for the high-dimensional model under control by Monte Carlo sampling with drastically reduced variance. We find that around 1000 particles are sufficient in practice. Our numerical experiments suggest that this method recovers the calibration speed from the corresponding 2-factor LSV model with deterministic rates defined in (2.2).

As a result of independent interest, we explain how to effectively deal with violation of the Feller condition for the Heston-type LSV Kolmogorov forward equation and numerically solve the PDE using a finite element method with a Backward Differentiation Formula (BDF) time-stepping scheme and an appropriate non-Dirichlet boundary condition. To the best of our knowledge, this represents a new approach which complements the literature on the use of ADI schemes [11, 39, 52] to handle the PDE calibration of an LSV model with deterministic rates.

For Heston type models, the CIR variance process can reach zero if the Feller condition is violated, as is often the case in FX markets (we refer to Table 6.5 in [11] for examples on a large range of currency pairs and maturities). As a consequence, the density is singular at the boundary $V=0$. In [45], the authors propose to reduce the problem by considering $\log\left(V_{t}/v_{0}\right)$, whereas [11] suggests to refine the mesh near $V=0$. While these methods alleviate the problem to some extent, we propose to use a different boundary condition as well as a change of variables which results in a bounded solution in a neighbourhood of $V=0$.

Moreover, a main advantage of the finite element method compared to ADI schemes, besides the greater flexibility in the mesh construction, is that the Dirac delta initial condition can be handled naturally in the weak formulation. This methodology yields an accurate calibration of the Heston-type LSV model with deterministic rates for a broad set of market data.

Finally, we provide empirical evidence that the inclusion of stochastic rates is important for the pricing of some specific exotic derivatives. In particular, in Section 5.4 we consider the pricing problem for a Target Accrual Redemption Note (TARN) and a no-touch option. We demonstrate that the impact of stochastic rates is comparable to the difference between pricing a 5-year no-touch option under a LV or LSV model. Other exotics with similar features, not considered here, are Accumulators and Power Reverse Dual-Currency notes (PRDC). Moreover, stochastic rates become necessary for any hybrid product which embeds the rates explicitly. Examples are spread options between an FX rate and the Libor rate.

The remainder of this paper is organised as follows. In Section 2, we specify the model and calibration framework and provide a necessary and sufficient condition for a perfect calibration to vanilla quotes. A rigorous proof emphasising the use of local times and possible moment explosions is given in Appendix A. In Section 3, we introduce the particle method used and detail how the control variates for both conditional expectations and standard expectations are constructed. In Section 4, we describe the calibration of the LSV model with deterministic rates using a carefully constructed finite element method. In Section 5, we present numerical results and show that a low number of particles suffices to provide a very good fit to market quotes, which demonstrates the computational efficiency of the method. The impact of stochastic rates for the pricing of a TARN and no-touch option is presented. Section 6 concludes with a brief discussion.

2 Model definition and calibration

We consider a domestic and a foreign market with stochastic short rates $r^{d}$ and $r^{f}$, and exchange rate $S$. The spot $S_{T}$ is associated with the currency pair ccy1ccy2 (following the notations in [11]) and denotes the amount of units of ccy2 (domestic currency) needed to buy one unit of ccy1 (foreign currency) at time $T$. We denote by $D^{d}$ and $D^{f}$ the domestic and foreign discount factors associated with their respective money market accounts,

[TABLE]

2.1 Models

We assume the existence of a filtered probability space ($\chi$, $\mathcal{F},\left\{\mathcal{F}_{t}\right\}_{t\geq 0},\mathbb{Q}^{d})$ with a domestic risk-neutral measure $\mathbb{Q}^{d}$. For future reference we also define a foreign risk-neutral measure $\mathbb{Q}^{f}$. Under $\mathbb{Q}^{d}$, $S$, $r^{d}$ and $r^{f}$ follow a system of SDEs

[TABLE]

where $V$ is the stochastic variance process and the four-dimensional standard Brownian motion $(W,W^{V},W^{d},W^{f})$ has the correlation structure

[TABLE]

with $\rho,\rho_{Sd},\rho_{Sf},\rho_{df}\in(-1,1)$, the other correlations being zero (and such that the correlation matrix is positive definite), and for given functions $\alpha:\mathbb{R}^{+}\times[0,T]\rightarrow\mathbb{R^{+}}$, $h^{d/f}:[0,T]\rightarrow\mathbb{R}$, and non-negative numbers $\kappa,\theta,\xi,\kappa_{d},\theta_{d},\xi_{d},\kappa_{f},\theta_{f},\xi_{f}$, as well as initial values $S_{0},g_{0}^{d},g_{0}^{f},V_{0}$.

Let the call option price under model (2.1) for a notional of one unit of $ccy1$, with strike $K>0$ and maturity $T>0$, be

[TABLE]

If the leverage function $\alpha\equiv 1$ in (2.1), we recover a Heston model with shifted CIR domestic and foreign short rates. We will refer to this model as Heston-2CIR++ model. As this model will only be used for intermediate calibration steps, we will make the additional simplification that the interest rate dynamics are independent of the dynamics of the spot FX rate and the variance process, for analytical tractability (see [3]).

We also define two simpler models which we will refer to in the remainder of the article. In both these models, rates are deterministic, $\bar{r}^{d}\left(t\right)=-\partial\ln P^{d}\left(0,t\right)/\partial t$ and $\bar{r}^{f}\left(t\right)=-\partial\ln P^{f}\left(0,t\right)/\partial t$, with $P^{d/f}\left(0,T\right)$ the market zero coupon bond prices for the domestic and foreign money market accounts, respectively.

We can then write the related 2-factor Heston-type LSV model with deterministic rates as

[TABLE]

for a given function $\alpha^{2D}:\,\mathbb{R}^{+}\times\left[0,T\right]\rightarrow\mathbb{R}^{+}$, and the pure Local Volatility (LV) model as

[TABLE]

with a given function $\sigma_{LV}:\,\mathbb{R}^{+}\times\left[0,T\right]\rightarrow\mathbb{R}^{+}$.

Note that while the volatility is “local”, i.e., a function of spot FX and time, the short rates are assumed to be a function of time only. We also note for future reference that under the pure LV model (2.3), call prices $C_{LV}$ satisfy the forward Dupire PDE (see [20])

[TABLE]

2.2 Calibration outline

The purpose of this paper is to calibrate $h^{d}$, $h^{f}$, $\kappa_{d}$, $\kappa_{f}$, $\theta_{d}$, $\theta_{f}$, $\xi_{d}$, $\xi_{f}$, $\kappa$, $\theta$, $\xi$, $\rho$ and especially $\alpha$ in (2.1). We will use calibration of (2.2) and (2.3) as “stepping stones”. More precisely, the full calibration process consists of the following steps, illustrated in Figure 2.1.

Calibration of Heston-2CIR++ parameters:

1. (a)

calibrate shifted CIR++ model for domestic and foreign short rates separately (Appendix D); 2. (b)

calibrate Heston-2CIR++ SV model assuming volatility, domestic and foreign short rates are all independent processes (Appendix F); 2. 2.

Calibrate local volatility assuming time-dependent domestic and foreign short rates (Appendix E); 3. 3.

Calibration of Heston-2CIR++ LSV model:

1. (a)

calibrate the leverage function of 2-factor Heston LSV model (2.2) assuming time-dependent domestic and foreign rates using local volatility from Step 2 (Section 4); 2. (b)

calibrate the leverage function of 4-factor LSV-2CIR++ LSV model (2.1) using Heston-2CIR++ parameters from Step 1, the local volatility from Step 2 and leverage function of 2-factor Heston LSV model as a control variate (Sections 3 and 5).

2.3 A necessary and sufficient condition for exact calibration

In the following, we give the main formula that links market call prices, via the Dupire local volatility, to prices under (2.1).

In [26], the following calibration condition is given111The (equivalent) context there is an equity with stochastic short rate and dividends.:

[TABLE]

where $\sigma_{LV}$ is a local volatility as in (2.3), and

[TABLE]

Assumption 1**.**

$\alpha$* is Lipschitz and uniformly bounded by $\alpha_{\max}$, $h^{d,f}$ are uniformly bounded and both the marginal density $\phi(\cdot,T)$ of $S_{T}$ in (2.1) and $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{T}^{d}V_{T}\,|\,S_{T}=\cdot\right]$ are continuous.*

We note that the continuity and positivity in $\mathbb{R}^{3}_{+}$ of the joint density $\psi(\cdot,\cdot,\cdot,T)$ of $\left(S_{T},V_{T},D_{T}^{d}\right)$ in (2.1) is sufficient for $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{T}^{d}V_{T}\,|\,S_{T}=\cdot\right]$ to be continuous.

We also define $\varphi=2+\sqrt{2}$, $\zeta=\xi\alpha_{\max}$, and

[TABLE]

which is a lower bound for the explosion time of $S_{t}^{2}$ (see [15]).

We prove the following theoretical results in Appendix A:

Proposition 2**.**

Under Assumption 1, the call price $C\left(K,T\right)$ under model (2.1) satisfies

[TABLE]

for any strike $K>0$ and maturity $T<T^{*}$, with $T^{*}$ given by (2.7).

Theorem 3**.**

Under Assumption 1, the call price $C\left(K,T\right)$ under model (2.1) matches the price $C^{LV}$ under the local volatility model (2.3) for any strike $K>0$ and maturity $T<T^{*}$ only if (2.5) holds for all $K,T>0$.

If (2.8) has a unique solution, then the condition (2.5) is also sufficient.

Uniqueness of the solution $C$ to the heat equation (2.8) is normally expected under sufficient regularity of the diffusion coefficient and under a growth condition.

The condition (2.5) expresses when a model of the form (2.1) with exogenously given $\alpha$ is consistent with market prices, which are expressed through the local volatility function $\sigma_{LV}$. We make no claim about the existence of such a model (see also Remark Remark below), and note that $\alpha$ enters (2.5) not only explicitly but also through the $\mathbb{Q}^{d}$-expectations. Existence of a calibrated model is linked to the existence of a solution to the McKean-Vlasov equation which results when inserting $\alpha$ defined endogenously by (2.5) in terms of $\sigma_{LV}$ and the model itself into (2.1). In [1], the existence of a short-time solution of the associated Fokker-Planck equation for the density of LSV processes of this type is shown under certain regularity assumptions. The upper bound on the time in [1] is needed to guarantee that the density stays strictly positive from an assumed strictly positive initial condition, and has no direct link to $T^{*}$ in this paper.

The ratio on the right-hand side of (2.5) accounts for the stochastic volatility; if there is no stochastic volatility (i.e. $V_{T}=1$), we recover the formula in [11]. The term $Q_{T}$ accounts for the stochastic rates and, if rates are deterministic, $Q_{T}=0$ and we recover the formula derived in [20],

[TABLE]

At time $T=0$, from $r_{0}^{d}-\bar{r}^{d}\left(0\right)=r_{0}^{f}-\bar{r}^{f}\left(0\right)=0$ we get $\alpha\left(K,0\right)=\sigma_{LV}\left(K,0\right)/\sqrt{v_{0}}.$

Theorem 3 provides technical conditions for the formula presented in [26], where a formal proof is given without specification of the rates processes. Here, we consider specifically an LSV–2CIR

++

model and derive the result rigorously. In Lemma 8 we provide a sufficient condition for the process

[TABLE]

to be a true martingale up to $T^{*}$, which is an important step in the proof of Theorem 3. On the one hand, $T^{*}$ is a lower bound for the explosion time of the second moment of the discounted spot process $D_{t}^{d}S_{t}$. On the other hand, [6] show that the moment explodes in finite time for the Heston model, a property that is inherited by our Heston-type LSV–2CIR

++

model (2.1) as well as the Heston-type LSV–2Hull–White model in [17]. Therefore, the formula may not hold for certain values of the model parameters and for large maturities $T$. However, in practice, $T^{*}$ is very large. For instance, from our calibration given in Section 5 we obtain $\kappa=1.4124$, $\xi=0.2988$, $\alpha_{\max}=1.40$, such that $T^{*}=28.6.$

Remark*.*

A numerical experiment in [26] raises the question of the existence of a calibrated 2-factor LSV model for large $\xi$ (there, $\xi\approx 350\%$ is used to match forward smiles). In this particular case and with the other model parameters kept the same, we find $T^{*}=0.20$, which indicates that moment explosions may occur sooner.

3 Fast calibration with a new control variate particle method

In this and the next section, we describe two of the main components of the calibration routine. We recall the calibration condition (2.5), which involves conditional expectations as well as standard expectations, which have to be estimated under model (2.1).

First, we describe the basic particle method for the estimation of these expectations. Then, we present the various control variates, building on intermediary calibration steps, which we use in order to reduce the computational cost of the calibration of $\alpha$ in the $4$-factor model (2.1).

Therefore, we require the prior calibration of the interest rate models in (2.1), the LV model (2.3), the Heston-2CIR++ model, and the LSV model (2.2). The calibration of the latter via a PDE is detailed in Section 4, while we refer to Appendices D, E, and F for the former three.

Equation (2.5) contains the local volatility, which can be obtained from derivatives of market prices from (2.4) by re-arranging it (into Dupire’s formula), and explicitly the second derivative of market prices with respect to strike. Different approximation approaches are used in practice, e.g., one writes the formulae in terms of the implied volatility, and uses a smooth parametrisation for the differentiation. Here, we first calibrate a parametrisation of the local volatility model with a fixed-point iteration as in [38, 46] and then use $\frac{\partial^{2}C_{LV}}{\partial K^{2}}$ obtained from the solution of the forward PDE (2.4) with a smoothing scheme (see Appendix E).

3.1 Calibration by particle method

A calibrated $\alpha$ is implicitly defined by (2.5), where the right-hand side depends on $\alpha$ in a non-linear way through the (conditional) expectations. Formal insertion of the calibration formula into the SDE (2.1) leads to a process where the diffusion coefficient depends on the distribution of the joint process $X_{t}=(S_{t},V_{t},r_{t}^{d},r_{t}^{f},D_{t}^{d})$. The process thus falls in the class of McKean-Vlasov processes [32].

The existence and uniqueness of the solution for this McKean-Vlasov SDE are not established theoretically, to the best of our knowledge. From an empirical perspective, in [26] and in Section 11.8 of [27] the authors encountered problems for very high values of $\xi$; see Remark Remark. In our case, for $\xi$ calibrated to market smiles ($\approx 30\%$) we are able to reach a high accuracy.

The particle method for processes of this type was introduced in [32] and is discussed in Chapter 2, Section 3 of [44]; it was applied to LSV model calibration in [26] and in Section 11.6 of [27].

We define $N$-sample path approximations of $X_{t}$ as $\left(X_{t}^{i,N}\right)_{i\in\llbracket 1,N\rrbracket}=(S_{t}^{i},V_{t}^{i},r_{t}^{d,i},r_{t}^{f,i},D_{t}^{d,i})_{i\in\llbracket 1,N\rrbracket}$ by the ($5\times N$)-dimensional SDE

[TABLE]

where $(W_{t}^{i},W_{t}^{d,i},W_{t}^{f,i},W_{t}^{V,i})$, ${i\in\llbracket 1,N\rrbracket}$ are $N$ i.i.d. copies of the four correlated Brownian motions, and $\hat{\alpha}_{N}$ is an estimator for $\alpha$ based on $\left(X_{t}^{i,N}\right)_{i\leq N}$,

[TABLE]

with

[TABLE]

where $\hat{p}_{N}$ is an estimator for $\mathbb{E}^{\mathbb{Q}^{d}}\!\!\left[D_{T}^{d}V_{T}\,|\,S_{T}=K\right]/\,\mathbb{E}^{\mathbb{Q}^{d}}\!\!\left[D_{T}^{d}\,|\,S_{T}=K\right]$, with $\delta_{N}$ a kernel function, and $Q_{T}^{i}$ is the $i$-th sample of $Q_{T}$ from (2.6) based on $X_{T}^{i,N}$.

The paths of the $5\times N$-dimensional process $(X_{t}^{i,N})_{i\leq N}$ are now entangled due to the dependence on $\hat{\alpha}_{N}$ in $(X_{t}^{i,N})_{i\leq N}$. The process can be seen as a system of $N$ interacting particles evolving in a $5$-dimensional space, where particle $i$ is defined by its position $X_{t}^{i,N}$. As in [26], we will therefore use the term “particle” instead of “path”. Because of the four driving factors, we will keep referring to this as a 4-factor model in spite of the extra state variable $D^{d}$.

A central ingredient for proving convergence of the particle method is the chaos propagation property (see Chapter 2, Section 3 of [44]), which is not proven for the present case.

3.2 Variance reduction for the Markovian projection

Our goal here is to reduce the variance of the estimator $\hat{p}_{N}$ from (3.3) to be able to use a minimal number of particles.

We assume that the 2-factor LSV model (2.2) is perfectly calibrated to market call prices, i.e. that (2.9) is satisfied. Then we will use

[TABLE]

which is an estimator for

[TABLE]

as a control variate for $\hat{p}_{N}$, and $p^{2D}$ will be computed using a PDE solver. The Kolmogorov forward equation for $p^{2D}$ is commonly used for the calibration of LSV models (see [11, 39, 52]), and we propose in Section 4 a new method which is tailored to the specific difficulties associated with density functions in Heston-style models.

We thus define a new estimator $p_{N}^{*}$ by

[TABLE]

The latter has an asymptotically diminishing bias if we assume the particle method to converge in distribution (and neglect the time stepping bias).

In order to get a good estimate for the optimal $\lambda$, we can rewrite the above estimator as

[TABLE]

with

[TABLE]

which mimics the standard Monte Carlo control variate form. We can think of $m_{i}$ and $m_{i}^{2D}$ roughly as samples of two random variables $m$ and $m^{2D}$ respectively (but note they are not independent, although for large $N$ the correlation is very low), and for the best variance reduction (see Section 4.1 in [23]), we take

[TABLE]

which we can estimate by

[TABLE]

We recall that the expected variance reduction factor is

[TABLE]

Hence, if the stochastic rates are not highly volatile, i.e. if $\xi_{f}$ and $\xi_{d}$ are small enough, we expect a very good variance reduction as the correlation between the particles generated by the 4-factor hybrid LSV model (2.1) and by the 2-factor LSV model (2.2) will be high. Our numerical tests performed on a model calibrated to recent EUR and USD market data exhibit a correlation between model (2.1) and (2.2) of $95\%$ up to $1.5$ years and $50\%$ around $5$ years. Additionally, our stress test in Subsection 5.3 suggests that even under high volatility regimes for the rate processes, i.e. when $\xi_{f}$ and $\xi_{d}$ are large, the variance reduction brought by this control variate is significant.

3.3 Variance reduction for standard expectations

Here, we discuss the variance reduction for the estimator $\widehat{Q}$ in (3.3) for $\mathbb{E}^{\mathbb{Q}^{d}}\left[Q_{T}\right]$, where we repeat

[TABLE]

from (2.6) for the convenience of the reader. This is an estimator for a standard expectation (in contrast to conditional expectations). We introduce control variates

[TABLE]

for $X_{1,T}$ defined in (3.8), and

[TABLE]

for $X_{2,T}$, and where $\widehat{r}^{f}$ is the foreign rate process without the quanto adjustment.222The last quantity is introduced because the expectation of ${r}^{f}$ is not analytically available. We know that if the model (2.1) is perfectly calibrated to call option prices,

[TABLE]

estimated from market data via a calibrated LV model. The following are also analytically available:

[TABLE]

We denote the Monte Carlo estimators of the corresponding $\mathbb{Q}^{d}$-expectations as $\widehat{X}_{1}$, $\widehat{X}_{2}$, $\widehat{Y}_{1}$, $\widehat{Y}_{2}$, $\widehat{Z}_{1}$, $\widehat{Z}_{2}$, respectively, using the same Brownian paths for $W,\;W^{V},\;W^{d},\;W^{f}$ in all estimators.

We can define a new Monte Carlo estimator ${Q}^{*}$ for $\mathbb{E}^{\mathbb{Q}^{d}}\left[Q_{T}\right]$ as

[TABLE]

with

[TABLE]

The weights $\lambda_{1}$, $\lambda_{2}$, $\eta_{1}$, $\eta_{2}$ above are chosen to minimize the variance of ${Q}^{*}$ (see [23]).

This approach is particularly useful for out-of-the-money options and digital options as the Monte Carlo estimator will exhibit higher variance in these settings.

A computational problem arises if there is no particle with $S_{T}>K$, which happens if $K$ is large and the total number of particles is relatively small (as will be the case with control variates), since the estimators of $\mathbb{V}\left[D_{T}^{d}\mathbf{1}_{S_{T}\geq K}\right]$ and $\mathbb{V}\left[D_{T}^{d}\left(S_{T}-K\right)^{+}\right]$ are then zero. In that case, we pick

[TABLE]

such that both control variates are of the same order of magnitude.

3.4 Implementation details

The leverage function $\alpha$ can in principle be computed for any $K$ and $T$ by the estimator (3.2). However, for computational purposes, we defined it in this way on a grid of points and interpolate it from there with cubic splines in spot and piecewise constant in time. We denote by $N_{T}$ the number of maturities. Then there are $N_{T}+1$ volatility “slices” in total such that we denote the $m$-th time slice $\alpha\left(\cdot,T_{m}\right)$, by $\alpha_{m}$, represented numerically as splines with $N_{S}$ nodes. While having $N_{S}$ too small will lead to accuracy problems, choosing it too large will make the surface rougher due to over-fitting. We find 25-30 points to provide a good trade-off between accuracy and smoothness. For a given $T_{m}$, the leverage function is thus defined on some interval $\left[S_{\min}^{m},\,S_{\max}^{m}\right]$ and is extrapolated constant outside these bounds. Because we need more grid points around the forward value and less around $S_{\min}^{m}$ and $S_{\max}^{m}$, we use a hyperbolic grid (with $\eta=0.05$, see Appendix C for more details) refined around the forward value

[TABLE]

where $\sigma_{F}\left(T_{m}\right)$ is the at-the-money forward market volatility for maturity $T_{m}$ (interpolated linearly in variance). Each of the grid values can be seen as a parameter and we denote them by $\left(\alpha_{m,j}\right)_{m\leq N_{T},\,j\leq N_{S}}$ with the associated spot grid values $\left(s_{m,j}\right)_{m\leq N_{T},\,j\leq N_{S}}$.

We now give the calibration algorithm. As previously, we denote the particle system at time $T$ for the model (2.1) by $\left(S_{T}^{i},V_{T}^{i},r_{T}^{d,i},r_{T}^{f,i},D_{T}^{d,i}\right)_{i\leq N}$. Similarly, we denote the 2-factor particle system at time $T$ for the model (2.2) by $\left(S_{T}^{2D,i},V_{T}^{2D,i}\right)_{i\leq N}$.

We work with an Gaussian kernel

[TABLE]

with a bandwidth given by a Silverman-type rule (see [42])

[TABLE]

where $\eta=1.5$ and $T_{\min}=0.25$ in our tests.

The step-by-step calibration is detailed in Algorithm 1.

Our empirical findings suggest that Quasi-Monte Carlo sampling of the random numbers does not provide significant accuracy gains. In order to speed up the computation of the sums involving kernel functions such as

[TABLE]

it is advised (see [26]) to sort the particle state vector by spot value and select only the relevant particles that fall inside an interval $\left[K-\Delta K,\,K+\Delta K\right]$, where we choose

[TABLE]

with $\epsilon=10^{-5}$.

4 Two-factor Heston-type LSV model calibration by PDE

Here, we describe the calibration of the 2-factor sub-model of (2.1) defined in (2.2) by solution of the forward PDE.333In this section only, we write $S$,$V$ and $\alpha$ in lieu of $S^{2D}$, $V^{2D}$ and $\alpha^{2D}$, respectively, to ease notation.

4.1 Transformation and weak formulation

The following is a small variation of the main result in [31], and the proof is therefore omitted. Note that a new non-Dirichlet boundary condition appears at $z=0$.

Theorem 4**.**

Define the region $\Omega=\mathbb{R}_{+}^{2}$ and assume that the density $\phi$ (under $\mathbb{Q}^{d}$) of the Markovian process $\left(S_{t},V_{t}\right)$ started at $\left(S_{0},v_{0}\right)$ at time 0 exists and is $C^{2,2,1}\left(\Omega\times\mathbb{R}_{+}\right)$. Then $\phi$ is the solution to the Kolmogorov forward equation

[TABLE]

Proof.

Similar to the proof of Lemma 4.1 and Theorem 4.1 in [31]. ∎

The marginal density function of the CIR process at $t$ is (see [14])

[TABLE]

with $c=\frac{2\kappa}{\left(1-e^{-\kappa t}\right)\xi^{2}}$ , $u=cv_{0}e^{-\kappa t}$ and $\beta=\frac{2\kappa\theta}{\xi^{2}}-1$, where $I_{\beta}$ is the modified Bessel function of the first kind of order $\beta$. We can write an asymptotic expression for small $z$ using the asymptotic formula for the modified Bessel function found in [2],

[TABLE]

such that $\phi_{v}\left(z,t\right)$ diverges for $z=0^{+}$ when $2\kappa\theta<\xi^{2}$. This agrees with the well-known density $\phi_{v}^{\infty}$ of the stationary distribution (see [14]) given by

[TABLE]

with $\omega=2\kappa/\xi^{2}$. This motivates a scaling of the density $\phi$ of the form $p=\phi z^{-\beta}$, and to solve a new PDE for $p$ which we define hereafter. By insertion in (4.1) we get the following.

Corollary 5**.**

For any $\beta\in\mathbb{R}$, $p=\phi z^{-\beta}$ satisfies the initial boundary value problem

[TABLE]

While this PDE is easier to handle numerically, one wants to work with the original density function $\phi$ for most of the applications. There are two main calculations one would like to achieve: the expected payoff $f\left(S_{T}\right)$ for a given function $f$; the Markovian projection $\mathbb{E}\left[V_{T}\,|\,S_{T}=K\right]$.

As $\phi$ is still intractable for small $z$ and computing $z^{\beta}p\left(x,z,t\right)$ is not numerically feasible, we perform an integration by parts (noticing $\lim_{z\rightarrow 0}z^{\beta+1}p=\lim_{z\rightarrow\infty}z^{\beta+1}p=0$ since $\beta+1>0$) to obtain, for deterministic rates,

[TABLE]

and

[TABLE]

4.2 Finite element method with a two-step BDF time scheme

We combine a finite element approximation in space with a Backward Differentiation Formula (BDF) scheme in time, since Crank-Nicholson time-stepping or ADI schemes can give rise to instabilities for Dirac initial data (see [36, 51]; we refer to [7] for a stability analysis and to [21] for some financial applications of BDF schemes).

Equation (4.4) can be written as

[TABLE]

Denote by $\Gamma_{R}=\left\{\left(x,z\right)\in\partial\Omega:\,z=0\right\}$ the subset of the boundary of $\Omega$ with Robin boundary condition. We derive a weak formulation in the usual way (see, e.g., [37]), i.e., we multiply the PDE by a test function $v\in H^{1}\left(\Omega\right)$, integrate over $\Omega$, using the divergence theorem and boundary conditions, to obtain the weak form of (4.2),

[TABLE]

with the bi-linear form

[TABLE]

where the last term contains the new boundary condition.

Let us define a uniform time mesh with $t_{m}=m\Delta_{t}$, $m\in\llbracket 0,M\rrbracket$. We denote $p_{m}=p\left(\cdot,\cdot,t_{m}\right)$, in which case the BDF scheme can be written as

[TABLE]

where the first time step is divided into two standard fully implicit time steps (see [22]). Then, for the first time step, using the Dirac delta initial condition,

[TABLE]

Note that we initially only assumed $v\in H^{1}\not\subset C$, and therefore the operation above with the Dirac delta is not defined for all such $v$. However, we will next use continuous basis functions. If $\left(S_{0},v_{0}\right)$ coincides with a mesh point, this is equivalent to solving a linear system where the right-hand side vector is $v_{0}^{-\beta}$ for the source point node and zero otherwise.

The PDE solution is approximated by a conforming finite element method with $\mathbb{P}_{2}$ elements, i.e., a polynomial of order two on a triangle cell. The computation of the finite element solutions was performed with the FEniCS library [4]. Each triangle is characterised by 6 local degrees of freedom (nodes) as displayed in Figure 4.2 (see [34] for details). We describe the mesh construction in detail in Appendix C. An example of a thus generated mesh with 30 spot steps and 30 variance steps is illustrated in Figure 4.3.

In order to see the improvement due to the transformed PDE (4.4) for $p$, we plot both $z^{-\beta}\phi$ and $p$ for a pure Heston model with

[TABLE]

which corresponds to a Feller ratio of $0.32$. We use 100 time steps as well as 30 spot steps and 30 variance steps. The solution for $p=z^{-\beta}\phi$ computed with no change of variables is presented in Figure 4.1, where we notice significant numerical instabilities. The tranformed PDE is solved for the same problem and $p$ is also plotted in Figure 4.1.

4.3 Calibration algorithm

The calibration of the 2-factor LSV model (2.2) is performed by finding the leverage function $\alpha$ defined in (2.9). We compute $\mathbb{E^{\mathbb{Q}}}\left[V_{T}\mid S_{T}=K\right]$ from (4.6) with the solution $p$ of the dampened PDE (4.4). Both integrals can be computed by double adaptive Clenshaw-Curtis quadrature rules to handle singularities properly when $T$ is small (see [24]).

Furthermore, for very small or very large values of $K$, both the numerator and denominator will be very small. So we define a smooth extrapolation rule by

[TABLE]

where we pick $\epsilon=10^{-14}$ in our numerical tests.

The calibration will be done forward in time. Denote by $\left(\Delta T\!\left(i\right)\right)_{i\leq N_{T}}$ the interval lengths between maturities and by $N_{T}$ the number of maturities.

The leverage function $\alpha$ is again defined by splines as detailed in Section 3.1. This approach allows us to compute $\mathbb{E^{\mathbb{Q}}}\left[V_{T}\,|\,S_{T}=K\right]$ only on the nodes, reducing the computational time considerably. In the calibration routine, we use forward constant interpolation of the leverage function between maturities to handle the non-linearity of the problem. We can then write the calibration procedure as in Algorithm 2.

In the test, we use the Heston parameters calibrated in Appendix F for the Heston-2CIR++ model, i.e.,

[TABLE]

where the Feller ratio is

[TABLE]

which violates the Feller condition.

The calibrated leverage function $\alpha$ is plotted in Figure 4.4. For the solution of the forward PDE between maturities, we use a BDF scheme with constant stepsize and find that $50$ time steps per year and a $80\times 80$ spot-variance finite element mesh give very accurate results.

5 Four-factor Heston-type LSV 2CIR++ model calibration

In this section, we give the results for the calibration of the main model (2.1) and test the efficiency of the algorithm.

We use vanilla options implied volatility data from Bloomberg from the $18/03/2016$ for the currency pair EURUSD, namely, $10\text{D-Put},25\text{D-Put},50\text{D},\,25\text{D-Call},\,10\text{D-Call}$, for the maturities444We skip the 7Y and 10Y quotes as they were not liquid enough. 3W, 1M, 2M, 3M, 6M, 1Y, 1Y6M, 2Y, 3Y, 5Y.

We use historical correlations estimated in [16] from weekly time series data from 2012–2014,

[TABLE]

We assume that both CIR

++

processes are calibrated under their own risk-neutral measure as in Appendix D, and that the Heston–2CIR

++

model is calibrated as in Appendix F for the Heston-2CIR++ model. The parameters are as follows

[TABLE]

In order to approximate the particle system (3.1), we use an extension of the QE-scheme from [5] to model (2.1). A full description of the time marching scheme is provided in Appendix B.

5.1 Calibration results and efficiency

For a first illustration of the model fit and the improvement through the control variates, we calibrate the 4-factor model with $800$ particles with and without control variates. The associated leverage function is plotted in Figure 4.4 (right).

We then plot, in Figure 5.1, the model implied volatility slices for 3M, 1Y, 2Y, and 5Y. The figure shows a significantly improved fit due to the control variates. We will analyse the accuracy and convergence in detail in Subsection 5.2.

As regards computational time, the control variate particle method with 20000 particles, which are enough for a reasonably converged solution (see Subsection 5.2), took approximately $10$% of the overall time spent in the PDE calibration of the LSV(2D) model. Most practitioners are familiar with the computational time required to calibrate a 2D LSV model by a forward PDE, a rough estimate being below one second for short- to mid-term expiries (5Y). The extra cost to calibrate the full 4D LSV–2CIR++ model is almost negligible with control variates, while the same accuracy without control variates requires more than 60 times the cost of the 2D PDE solver.

5.2 Variance and error reduction

In this subsection, we compare the results with control variates (CV) to those with the plain particle method (No-CV) as a function of the number of particles. The error measure we use is the absolute error in volatility (in $\%$ units). For instance, a maximum error (taken over all quoted deltas and maturities) of $0.03\%$ for a $20\%$ market volatility signifies that the calibrated volatility can be $20.00\%\pm 0.03\%$ in the worst case scenario.

We use the calibration routine described in Subsection 3.1 with and without control variates for 160, 800, 4 000, 20 000, 100 000, 500 000 and 2 500 000 particles.

The results are presented in Figure 5.2. On the basis that higher correlations between spot-rates and rate-rate will make our control variates more efficient, we also display in Figure 5.2 the calibration errors with no correlations between spot-rates and rate-rate as a presumed worst-case.

We infer from the data in Figure 5.2 that the use of control variates greatly improves the general calibration routine. Convergence in both the maximum error ($0.032\%$) and in the average error ($0.012\%$) is reached (i.e., the error from here on is dominated by other sources, such as the time discretisation error) for 4 000 particles when using control variates, while the plain particle method (without variance reduction) only reaches the same accuracy with 2 500 000 particles. We infer that the control variates give a $625$-fold speed-up.

Moreover, from careful data analysis of Figure 5.2, we find a convergence rate of $0.3$ in the number of particles for the plain particle method, for both the average and maximum error. The addition of the control variates preserves the convergence rate but reduces the absolute size of the error significantly.

We note that the error can be further reduced by increasing the number of time steps per year. We used $250$ time steps (such that, on average, there is one time step per open day) as it already yields a very accurate calibration at a reasonable computational cost.

Part of the error reduction is due to the conditional control variate described in Subsection 3.2 which can provide very good results for short-term horizons (the other part being due to the control variates for standard expectations). In order to analyse this further, we plot the variance reduction factor from (3.7), estimated with 500 000 particles, as a function of time in Figure 5.3. We estimate a trend line $1+\frac{C}{T^{1.6}}$, for a given constant $C$, and thus a very good variance reduction for short-term options. For longer terms, this still yields a good variance reduction factor of $5.2$ for the $3Y$ maturity and $2.34$ for $5Y$.

In addition, the other two control variates presented in Subsection 3.3 help control the short- to long-term behaviour as well. We plot both variance reduction factors for $X_{1}^{*}$ and $X_{2}^{*}$ in Figure 5.3. They seem to reach a steady state for maturities around $1.6$ and $6.5$, respectively. We note that we displayed here average variance reduction values over all strikes, while $X_{1}^{*}$ was found to provide significant variance reduction for small strikes.

5.3 Stress scenario

The last calibration result we present is for an extremely stressed set of rates parameters. This will serve as a robustness test of the control-variate particle method for very volatile short rate processes. In order to perform the test, we multiply the calibrated $\xi_{d}$ and $\xi_{f}$ by $20$ as well as divide $\theta_{d}$ and $\theta_{f}$ by $20$. This leads to strongly violated Feller conditions and very high volatilities for the two CIR

++

rate processes. We also set $\rho_{Sd}=0,\rho_{Sf}=0,\rho_{df}=0$ as this has shown to be more challenging for our method (as it makes the control variates less effective). The stress scenario parameter values are

[TABLE]

We emphasise that in practice, the volatility parameters $\xi_{d/f}$ are rarely above $0.06$. We refer the reader to [9] for more details. A calibration summary for the average error in absolute volatility is displayed in Table 1. In this stress scenario, the calibration via control-variate particle method reaches an error of $0.0165\%$ for $20\,000$ particles, whereas $2\,500\,000$ particles are required with the plain particle method (without control variates) to reach the same accuracy. Hence, the control-variate particle method shows a consistent improvement over the plain particle method even under stress scenarios. The conditional control variate for $p_{N}^{*}$ yields a variance reduction factor of almost $2$ for the last maturity ($5\text{Y}$), similar to the correlated case.

5.4 Impact of stochastic rates

In this subsection, we discuss the pricing of more exotic products, namely a no-touch option and a target accrual redemption note (TARN); see Chapter 8 in [11] and Section 2.2 in [53] respectively for a discussion of these products. Specifically, we assess the impact of stochastic rates on products embedding knock-out features with mid- to long-term expiries.

No-touches

The foreign no-touch up option pays one EUR at maturity if the exchange rate has not breached an upper barrier during the product lifespan. The payout under an arbitrage-free model with foreign risk-neutral measure $\mathbb{Q}^{f}$ , for a premium expressed in foreign units (EUR), a notional $N_{\text{EUR}}$ in foreign units (EUR) and for a maturity $T$ is

[TABLE]

where $M_{t}=\sup_{0\leq u\leq t}S_{u}$ is the running-maximum of the spot $S$ and $B_{1}$ is the upper barrier. For the tests, we pick

[TABLE]

TARNs

The specification of the TARN is here as follows: the buyer receives the forward value $K-S_{t_{i}}$ at fixing date $t_{i}$ if $K>S_{t_{i}}$; the buyer has to pay $S_{t_{i}}-K$ if $K<S_{t_{i}}$. At each fixing date, if the amount received is positive, the accrued value is increased by the paid amount. If at some point in the deal life-cycle, the accrued amount breaches a target $H_{target}$, the deal is terminated early. To protect the buyer, an additional knock-out barrier redeems the deal early if the spot fixes at or above an upper barrier $B_{2}$. The payout under an arbitrage-free model with domestic risk-neutral measure $\mathbb{Q}^{d}$ , for a premium expressed in domestic units (USD) and a notional $N_{EUR}$ in foreign units (EUR), for a maturity $T$ is

[TABLE]

where $\tau$ is the (early) redemption date that is either triggered by the accrual breaching the target $H_{target}$ or the spot $S$ breaching the upper barrier $B_{2}$. In this test, we pick monthly fixings and

[TABLE]

We display in Table 2 the prices computed with 1 048 575 quasi-Monte Carlo paths (with Sobol sequences and Brownian bridge construction from [10]) and 365 time steps per year for the foreign no-touch and TARN contracts under the LV, LSV and LSV–2CIR++ models. The running maximum is sampled with the Brownian bridge technique, as described in Chapter 6 of [23].

We infer from the data in Table 2 a relative difference of $1.35\%$ in price for the foreign no-touch and $-16.07\%$ for the TARN between the LSV and LSV–2CIR

++

models, which confirms that adding stochastic rates has a significant impact even for a 5Y contract.

Interestingly, in this setup, the price of the foreign no-touch option under the LSV–2CIR

++

lies in between the LV and LSV models. It is common to observe a market price that lies between the LV model price and a Heston-type LSV model price, when the Heston parameters have been calibrated to the same vanilla options. Practitioners therefore introduce a mixing factor to control the amount of stochastic volatility of the LSV model and to manually match the no-touch options quotes (see [11] for details). The introduction of stochastic rates seems to achieve a similar behaviour, at least in this example.

6 Conclusion

In this paper, we have provided a new and numerically effective method to calibrate a 4-factor LSV model to vanilla options. In our numerical tests with market data, we managed to achieve an approximate $625$-fold speed-up for the calibration using control variates, as compared to the plain particle method. We have shown that a high accuracy can be obtained with as few as 4 000 particles (with a maximum error in absolute volatility of $0.03\%$), and we were able to get a good fit with only 800 particles (with a maximum error in absolute volatility of $0.05\%$).

Using the calibrated leverage function from this paper, we showed that the addition of stochastic rates has a significant impact on structured products, even more so when barrier features and coupon detachments are combined for longer-dated contracts. Stochastic rates become necessary in the modelling if one wishes to price hybrid products where the rates appear explicitly (for instance, a spread option on the FX performance and the Libor rate). One could use a second factor in the CIR

++

processes to improve the fit to caps and use the method presented in this article to calibrate the leverage function.

Appendix A Proof of Proposition 2

and Theorem 3

We first state a necessary auxiliary result, which is an adaptation of Tanaka’s formula in Chapter 4 of [40], where the integrand $D$ in the local time integral is $1$.

Proposition 6**.**

On a filtered probability space ($\chi$, $\mathcal{F},\left\{\mathcal{F}_{t}\right\}_{t\geq 0},\mathbb{Q}^{d})$, let $D$ and $X$ be two $\mathcal{F}_{t}$-adapted continuous semi-martingales, with $D$ positive and integrable, $\left(l_{t}^{a}\right)_{t\geq 0}$ the local time of $X$ at level a and, for all $n>0$,

[TABLE]

Then, for any $T>0$ and $0<t\leq T$, $\int_{0}^{t}D_{s}\delta_{n}^{a}\left(X_{s}\right)\,d\langle X\rangle_{s}$ converges almost surely, and uniformly in time, to $\int_{0}^{t}D_{s}\,dl_{s}^{a}$.

Proof.

We closely follow Section 45 of [40], but use a more concrete expression of the regularisation function. The local time at level $a$ is defined as a continuous adapted increasing process such that

[TABLE]

with $\operatorname{sgn}\left(x\right)=-1$ for $x\leq 0$ and $\operatorname{sgn}\left(x\right)=1$ for $x>0$. We define a sequence of functions $\left(f_{n}\right)_{n\geq 0}$ for all $x\in\mathbb{R}$, as in Chapter 4 of [35]

[TABLE]

Hence, for all $n>\text{0}$, $\frac{1}{2}f_{n}^{{}^{\prime\prime}}=\delta_{n}^{a}$ a.e. We recall from the proof of Tanaka’s formula in [40] that $\int_{0}^{t}\delta_{n}^{a}\left(X_{s}\right)\,d\left\langle X\right\rangle_{s}$ converges almost surely to $l_{t}^{a}$ (uniformly in $t$). We note that the sequence $f_{n}$ converges uniformly to $x\rightarrow\left|x-a\right|$ and $f_{n}^{{}^{\prime}}$ converges point-wise to $\operatorname{sgn}\left(x-a\right)$. By the Itô-Doeblin formula we can write

[TABLE]

We denote

[TABLE]

and, since $f_{n}^{{}^{\prime\prime}}\left(x\right)=0$ for any $x$ such that $\left|x-a\right|\geq\frac{1}{n}$, we have

[TABLE]

Also, from the definition of $f_{n}$, for all $x\in\mathbb{R}$,

[TABLE]

and then for any given $x$, and $n>0$, both $D_{t}\left|\operatorname{sgn}\left(x-a\right)-f_{n}^{{}^{\prime}}\left(x\right)\right|$ and $D_{t}\left|\left|x-a\right|-f_{n}\left(x\right)\right|$ are smaller than $D_{t}$ which is integrable. Let $X_{t}=X_{0}+M_{t}+A_{t}$ be the canonical decomposition of $X$ and $D_{t}=D_{0}+N_{t}+R_{t}$ the canonical decomposition of $D$. Localisation allows us to reduce the problem to the case where $M$ and $N$ are bounded and $A$ and $R$ are of bounded variation. Then,

[TABLE]

for which the right-hand-side goes to zero when $n$ goes to infinity. By Doob’s $L^{2}$ martingale inequality, we can write

[TABLE]

to conclude that

[TABLE]

in $L^{2}$ and in probability. We may then assume that (A.4) also holds almost surely (since we could work with a subsequence for which the statement is true instead). Similarly, we have

[TABLE]

Also,

[TABLE]

and

[TABLE]

The monotone-convergence theorem allows us to conclude that

[TABLE]

in $L^{1}$, in probability and almost surely (on passing to a sub-sequence on $n$ if necessary). Hence,

[TABLE]

almost surely and uniformly in time. Similarly,

[TABLE]

and we get

[TABLE]

almost surely and uniformly in time. Additionally, we can write

[TABLE]

From the Kunita-Watanabe inequality,

[TABLE]

Since $\left\langle N\right\rangle_{s}$ and $\left\langle M\right\rangle_{s}$ are increasing processes of finite variation, we proceed as in (A.5) and conclude

[TABLE]

Hence, from (A), $C^{n}$ converges to a limit $\zeta$ almost surely (uniformly in time $t$). Applying integration by parts to the Tanaka formula (A.2), we can write

[TABLE]

where Tanaka’s formula also allows us to write

[TABLE]

and conclude that

[TABLE]

∎

We now state and prove two lemmas necessary for the derivation of Theorem 3. Lemma 7 provides a link between the local time of a process and its density function.

Lemma 7**.**

Given a filtered probability space ($\chi$, $\mathcal{F},\left\{\mathcal{F}_{t}\right\}_{t\geq 0},\mathbb{Q}^{d})$, let $W$ be a standard Brownian motion and $\mu$, $Y$ two $\mathcal{F}_{t}$-adapted processes with $Y$ continuous, with finite second moment and

[TABLE]

Consider a continuous Itô process $X$ given by

[TABLE]

whose marginal density function $\phi(\cdot,t)$ and $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{t}Y_{t}^{2}\,|\,X_{t}=\cdot\right]$ are assumed to be continuous. Further denote by $\left(l_{t}^{a}\right)_{t\geq 0}$ the local time of $X$ at level $a$. Then, for any continuous, integrable and positive $\mathcal{F}_{t}$-adapted semi-martingale $D$ and any $a\in\mathbb{R}$,

[TABLE]

Proof.

From Proposition 6, we know that the following holds almost surely:

[TABLE]

This implies convergence in distribution, so that we can write

[TABLE]

and, by the stochastic Fubini theorem, we get

[TABLE]

We denote $\gamma\left(x,u\right)=\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{u}Y_{u}^{2}\,|\,X_{u}=x\right]$, such that

[TABLE]

where we have used Fubini’s Theorem in the second line. By the continuity assumptions on $\phi$ and $\gamma$, we deduce that

[TABLE]

∎

Lemma 8**.**

Given the set-up of Theorem 3,

[TABLE]

is a true martingale up to $T^{*}$ given by (2.7).

Proof.

Since $\alpha$ and $\mathbf{1}_{S_{u}\geq K}$ are bounded, the process

[TABLE]

is a true martingale if

[TABLE]

On the one hand, since $t<T^{*}$, from Proposition 3.13 in [15], we can find $\omega>2$ such that

[TABLE]

On the other hand, from Theorem 3.1 in [29],

[TABLE]

Using Hölder’s inequality with the pair $\left(\frac{\omega}{2},\frac{\omega}{\omega-2}\right)$,

[TABLE]

Finally, using the Fubini theorem, $\mathbb{E}^{\mathbb{Q}^{d}}\left[\int_{0}^{t}\left(D_{u}^{d}S_{u}\right)^{2}V_{u}\,du\right]<\infty$ and hence $M$ is a true martingale of zero expectation. ∎

Combining Lemmas 7 and 8, we can derive Proposition 2.

Proof of Proposition 2.

Let $K\in\mathbb{R}^{+}$, $0<t<T^{*}$ and $H_{t}=\left(S_{t}-K\right)^{+}$. The Trotter-Meyer theorem [40] gives

[TABLE]

which we can write in differential form as

[TABLE]

Also,

[TABLE]

Hence, by applying Lemma 7 with $X_{t}=S_{t}$, $D_{t}=D_{t}^{d}$ and $Y_{t}=\alpha\left(S_{t},t\right)S_{t}\sqrt{V_{t}}$, we can write

[TABLE]

where $\phi$ is the marginal density function of $S$ at time $t$. Furthermore, one can define $\bar{\phi}_{n}$ as

[TABLE]

with $\delta_{n}^{K}$ defined as in (A.1) by

[TABLE]

and by a similar reasoning to that of Lemma 7 we get

[TABLE]

Since

[TABLE]

we write

[TABLE]

Combining (A.7) and (A.8) allows to write

[TABLE]

Hence, integrating (A.6),

[TABLE]

Furthermore, on a fixed time interval $\left[0,T^{*}\right]$, $D_{t}^{d}$ is uniformly bounded by $\exp({T^{*}\max_{u\in[0,T^{*}]}\left|h^{d}\left(u\right)\right|})$. Then, from Lemma 8 we know that $\int_{0}^{t}\mathbf{1}_{S_{u}\geq K}D_{u}^{d}\alpha\left(S_{u},u\right)S_{u}\sqrt{V_{u}}\,dW_{u}$ is a true martingale of zero expectation.

We write (A) at time $T$, differentiate with respect to $T$ and, upon noticing that $\mathbf{1}_{S_{T}\geq K}S_{T}=\left(S_{T}-K\right)^{+}+\mathbf{1}_{S_{T}\geq K}K$, we get (2.8). ∎

We are now ready to give the proof of Theorem 3.

Proof of Theorem 3.

First, we want to ensure that (2.5) is a necessary condition for

[TABLE]

Hence, by subtracting the Dupire PDE (2.4) from (2.8), we obtain

[TABLE]

so

[TABLE]

where

[TABLE]

It remains to show that we can replace $\overline{Q}_{T}$ by $Q_{T}$ in $\bar{q}$ . First, $D_{T}^{d}\left(S_{T}-K\right)^{+}$ is weakly differentiable with respect to $K$ with $\frac{\partial D_{T}^{d}\left(S_{T}-K\right)^{+}}{\partial K}=D_{T}^{d}\mathbf{1}_{S_{T}\geq K}$, which is bounded by the integrable process $D_{T}^{d}$. We can interchange differentiation and expectation to get $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{T}^{d}\mathbf{1}_{S_{T}\geq K}\right]=-\frac{\partial C}{\partial K}$. Since the models agree, $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{T}^{d}\mathbf{1}_{S_{T}\geq K}\right]=-\frac{\partial C_{LV}}{\partial K}$ and $\mathbb{E}^{\mathbb{Q}^{d}}\left[D_{T}^{d}\left(S_{T}-K\right)^{+}\right]=C_{LV}$, and (2.5) holds.

By re-tracing the steps in reverse order, one sees that (2.5) is also a sufficient condition for (A.10) provided the solution to (2.8) is unique. ∎

Appendix B Monte Carlo QE-scheme

The Quadratic-Exponential (QE) scheme [5] that is used to discretise the square-root process, employs moment-matching techniques and can significantly reduce the Monte Carlo discretisation error. While the full truncation Euler scheme and the QE scheme have shown to perform well in our tests, we experienced a faster convergence in time for the QE scheme when the Feller condition is broken. Hence, we choose the QE scheme for the variance process and the full truncation Euler for both stochastic rates, as the computational cost will be smaller. We briefly write a generalised QE scheme based on the original scheme from [5] to incorporate a leverage function and stochastic rates in the discretisation.

We follow our time interpolation rule for the calibration of $\alpha$ and interpolate forward-flat in time. We assume for simplicity that each Monte Carlo time step belongs to the $\alpha$ time grid. We can write

[TABLE]

and hence

[TABLE]

and

[TABLE]

where $W_{t}^{S}$ is a Brownian motion independent of $W_{t}^{V}$. Therefore,

[TABLE]

We approximate $\int_{t}^{t+\Delta t}V_{u}\,du$ by $\left(\frac{V_{t+\Delta t}+V_{t}}{2}\right)\Delta t$ , and note that conditional on $V_{t}$ and $\int_{t}^{t+\Delta t}V_{u}\,du$, since $W_{u}^{V}$ and $W_{u}^{S}$ are independent, the Itô integral $\int_{t}^{t+\Delta t}\sqrt{V_{u}}\,dW_{u}^{S}$ is normally distributed with mean zero and variance $\int_{t}^{t+\Delta t}V_{u}\,du$. We write the full scheme below

[TABLE]

with

[TABLE]

Let the Cholesky decomposition of the correlation matrix555If the correlation matrix is not positive definite, one can rely on spectral decomposition instead; see 2.3 in [23] for details.

[TABLE]

be $\mathbf{LL}^{T}$. $Y_{d}$ and $Y_{f}$ are defined as

[TABLE]

where $Z,\,Z_{v},\,Z_{d},\,Z_{f}$ are independent draws from a standard normal distribution and $U$ is a draw from a uniform distribution.

Appendix C Finite element mesh construction

In order to refine the mesh in the most relevant area, we use an exponential mesh on the variance axis and a hyperbolic mesh (see [50]) in the spot direction. This makes the mesh finer around $z=0$ and $x=S_{0}$. In order to build our mesh, we first define the grids in spot \left(x_{i}\right)_{i\in\left\llbracket 0,N_{S}\right\rrbracket} and variance \left(z_{j}\right)_{j\in\left\llbracket 0,N_{V}\right\rrbracket} separately. Additionally, to solve the PDE numerically, we need to truncate at the boundary and use $\bar{\Omega}=\left\{\left(x,z\right)\in\left[0,S_{\max}\right]\times\left[0,V_{\max}\right]\right\}$ on a time interval $\left[0,T\right]$ . We choose

[TABLE]

and recall that the stationary distribution of the CIR process is a gamma distribution of density $\phi_{v}^{\infty}$ as defined in (4.3). We compute $V_{\max}$ with the inverse cumulative density function such that

[TABLE]

We write

[TABLE]

with

[TABLE]

where $\eta$ is defined according to our numerical experiments and $g_{h}$ is the quadratic polynomial that passes through the points $\left(0,0\right)$, $\left(1,1\right)$,

[TABLE]

[TABLE]

where $\lambda$ is defined according to our numerical experiments and $g_{h}$ is the quadratic polynomial that passes through the points $\left(0,0\right)$, $\left(1,1\right)$,

[TABLE]

The latter intermediate step makes sure that both $S_{0}$ and $v_{0}$ are vertices of their respective grids. The construction of the finite element triangular mesh can be achieved by creating a vertex at each point $\left(x_{i},z_{j}\right)$ and defining two triangular cells (upper left and lower right) in each rectangle.

Appendix D Shifted CIR model and calibration

The domestic and foreign short interest rates are modeled by the shifted CIR (CIR

++

) process [9]. On the one hand, this model preserves the analytical tractability of the CIR model for bonds, caps and other basic interest rate products. On the other hand, it is flexible enough to fit the initial term structure of interest rates exactly. For $i\in\{d,f\}$, the short rate dynamics under their respective spot measures, i.e., $\mathbb{Q}^{d}$ – domestic and $\mathbb{Q}^{f}$ – foreign, are given by

[TABLE]

where $B^{d}$ and $B^{f}$ are Brownian motions under $\mathbb{Q}^{d}$ and $\mathbb{Q}^{f}$, respectively. The mean-reversion parameters $\kappa_{i}$, the long-term mean parameters $\theta_{i}$ and the volatility parameters $\xi_{i}$ are the same as in (2.1). The calibration of the short rate model (D.1) follows the same approach for both the domestic and the foreign interest rate. For simplicity, we drop the subscripts and superscripts “$d$” and “$f$” in the remainder of the subsection and define the vector of parameters $\beta_{1}=(g_{0},\kappa,\theta,\xi)$. According to Brigo and Mercurio [8], an exact fit to the initial term structure of interest rates is equivalent to $h(t)=\varphi^{\text{CIR}}(t;\beta_{1})$ for all $t\in[0,T]$, where

[TABLE]

$\nu=\sqrt{\kappa^{2}+2\xi^{2}}$ and $\bar{r}(0,t)$ is the market instantaneous forward rate at time [math] for a maturity $t$, i.e.,

[TABLE]

where $P(0,t)$ is the market zero coupon bond price at time [math] for a maturity $t$. The value of the zero coupon bond is given by

[TABLE]

where $\Delta(0,t)$ is the year fraction from [math] to time $t$ and $R(0,t)$ is the current (simply-compounded) deposit rate with maturity date $t$ which is quoted in the market. As an aside, note that the standard day count convention for USD and EUR is Actual $360$.

The detailed calibration procedure for both domestic and foreign rate processes can be found in Appendix D.1. The calibration results are displayed in Table 3.

D.1 Shifted CIR model calibration

In order to estimate the zero coupon curve (also known as the term structure of interest rates or the yield curve), we assume that the instantaneous forward rate is piecewise-flat. Consider the time nodes $t_{0}\!=\!0,t_{1},\ldots,t_{n}$ and the set of estimated instantaneous forward rates $f_{1},f_{2},\ldots,f_{n}$ from which the curve is constructed, and define

[TABLE]

Using (D.3) – (D.5) and solving the resulting linear system of equations, we get

[TABLE]

The continuously-compounded spot rate, i.e., the constant rate at which the value of a pure discount bond must grow to yield one unit of currency at maturity, is defined as

[TABLE]

Using (D.5) – (D.7), we deduce that

[TABLE]

whenever $t_{i-1}\leq t<t_{i}$. In Figure D.1, we plot the USD and EUR zero coupon curves $t\mapsto R_{0}(0,t)$, $t>0$, estimated from the quoted deposit rates from March 18, 2016, together with the flat-forward instantaneous forward rates.

A choice of the shift function $h$ as in (D.2) results in an exact fit to the initial term structure of interest rates independent of the value of the parameter vector $\beta_{1}$.

Next, we determine $\beta_{1}$ by calibrating the CIR

++

model to the current term structure of volatilities, in particular, by fitting at-the-money (ATM) cap volatilities. We consider caps with integer maturities ranging from $1$ to $10$ years for both currencies, with an additional $18$ month cap for EUR. For USD, all caps have quarterly frequency, whereas for EUR the $1$ year and $18$ month caps have quarterly frequency and the $2$ to $10$ year caps have semi-annual frequency. A cap is a set of spanning caplets with a common strike so the value of the cap is simply the sum of the values of its caplets. It is market standard to price caplets with the Black formula, in which case the fair value of the cap at time [math] with rate (strike) $K$, reset times $T_{a},T_{a+1},\ldots,T_{b-1}$ and payment times $T_{a+1},\ldots,T_{b-1},T_{b}$ is:

[TABLE]

where $F(0,T,S)$ is the simply-compounded forward rate at time [math] for the expiry $T$ and maturity $S$ defined as

[TABLE]

and the Black volatility $\sigma_{a,b}$ corresponding to a strike $K$ is retrieved from market quotes. Denoting by $\phi_{0}$ and $\Phi_{0}$ the standard normal probability density function (PDF) and cumulative distribution function (CDF), respectively, Black’s formula is:

[TABLE]

However, Black’s formula cannot cope with negative forward rates $F$ or strikes $K$, in which case we switch to Bachelier’s (normal) formula in (D.9):

[TABLE]

The data in Figure D.1(b) suggest that the instantaneous forward rate for EUR takes negative values. Therefore, we use Black cap volatility quotes for USD and Normal cap volatility quotes for EUR. The market prices of at-the-money caps are computed by inserting the forward swap rate

[TABLE]

as strike and the quoted cap volatility as $\sigma_{a,b}$ in (D.9), using either Black’s or Bachelier’s formula.

Fitting the CIR

++

model to cap volatilities means finding the value of $\beta_{1}$ for which the model cap prices, which are available in closed-form [9], best match the market cap prices. The calibration is performed by minimising the sum of the squared differences between model- and market-implied cap volatilities:

[TABLE]

where $\sigma^{\text{CIR}}$ and $\sigma^{\text{M}}$ stand for the model- and the market-implied cap volatilities, respectively, and $T_{1},\ldots,T_{n}$ are the cap maturities. Model-implied cap volatilities are obtained by pricing market caps with the CIR

++

model and then inverting the formula (D.9) in order to retrieve the implied volatility associated with each maturity. We choose to calibrate the model to cap volatilities since they are of similar magnitude, unlike cap prices which can differ by a few orders of magnitude. The calibration results are displayed in Table 3.

On the one hand, (D.14) is a highly nonlinear and non-convex optimisation problem, and the objective function may have multiple local minima. On the other hand, global optimisation algorithms require a very high computation time and do not scale well with complexity, as opposed to local optimisation methods. A fast calibration is important in practice since option pricing models may need to be re-calibrated several times within a short time span. Therefore, we used a nonlinear least-squares solver, in particular the trust-region-reflective algorithm [13], for the calibration and a global optimisation method, in particular a genetic algorithm [43], for verification purposes only.

Figure D.2 shows the fitting capability of the CIR

++

model, and the implied cap volatility curve is compared to the market curve for each currency. Taking into account that the model has only $4$ parameters to fit between $10$ and $11$ data points, we conclude that the CIR

++

model provides a fairly reasonable fit to the term structure of cap volatilities $T_{i}\mapsto\sigma^{\text{M}}(T_{i})$, $1\leq i\leq n$.

Appendix E Local volatility calibration algorithm

E.1 Calibration with Dupire PDE

The calibration routine for a pure local volatility model is run with a standard algorithm forward in maturity. We recall that model (2.3) is written as

[TABLE]

and we want to find the function $\sigma_{LV}$ for which the call prices under the local volatility model match the quoted market prices exactly. This is crucial as both $\sigma_{LV}$ and $\frac{\partial^{2}C_{LV}}{\partial K^{2}}$ appear in the leverage function formula (2.5). The forward Dupire PDE (2.4),

[TABLE]

provides an efficient way to calibrate and, eventually, regularise the problem. Denote by $\Phi$ the map from the local volatility function to the model implied volatility function $\Sigma_{Model}$. Furthermore, the PDE solution for a guess $\sigma_{S}$ of the local volatility gives call prices for the whole set of strikes and maturities. Inverting the Black formula allows to retrieve the model implied volatilities $\Sigma$. Hence, as proposed in [46], we can use the forward Dupire PDE (E.1) combined with an efficient implied volatility inverter [30] as the mapping function $\Phi$. A very useful property of this PDE is that it can be solved forward in maturity. Let a set of maturities quoted on the market be $\left(T_{1},...,T_{N_{Mat}}\right)$ and a set of $M_{i}$ strikes for a given maturity $T_{i}$ be $\left(K_{T_{i},1},...,K_{T_{i},M_{i}}\right)$. It is possible to solve the PDE on $\left[0,T_{1}\right]$, then on $\left[T_{1},T_{2}\right]$ and so forth. The full calibration algorithm is presented for completeness in Appendix E.

E.2 Computation of the target volatility surface

For the calibration routine, we will compute the solution of the PDE (E.1) by a finite difference method. The spot grid is defined on $\left[0,S_{\max}\right]$, where $S_{\max}=S_{0}e^{\frac{6}{2}\sigma_{ATM}\left(\frac{T_{\max}}{2}\right)\sqrt{\frac{T_{\max}}{2}}}$. In order to speed up the calibration routine, we prefer not to use too many spot steps and time steps (150 steps in space and 20 time steps per year). Hence, the scheme will not have converged to the solution of the PDE at this point. In order to tackle this problem and still benefit from a good speed-up, we will compute a “target volatility surface”: instead of calibrating the market volatility surface, we will calibrate a volatility surface that takes into account the discretisation error of the numerical PDE solution. Industry practitioners like Murex use this approach [33]. The algorithm to build the target surface is explained below.

E.3 Calibration by fixed-point algorithm and forward induction

The local volatility function is defined on a grid of points interpolated with cubic splines in spot and backward flat in time. In the FX case, where there are 5 quoted strikes per maturity (10 maturities), the local volatility is defined on a grid of $50$ points. Each one of the points $\sigma_{LV}^{i,j}=\sigma_{LV}\left(K_{T_{i},j},T_{i}\right)$, with $i\in\llbracket 1,10\rrbracket$ and $j\in\llbracket 1,5\rrbracket$, can be seen as a parameter of the local volatility surface. For a given maturity $T_{i}$, the local volatility is defined on the interval $\left[K_{T_{i},1},K_{T_{i},5}\right]$ and is extrapolated flat outside those bounds.

In order to define a first guess for the calibration routine, we use a smoothed bi-variate cubic spline following the algorithm in [19] to interpolate in strike and maturity the call prices on the market. This allows us to use the Dupire formula to define a first guess for the first maturity $T=T_{1}$. After the calibration of the first maturity pillar $T_{1}$, the first guess for the next pillar is the current maturity local volatility. This approach has shown the best stability and speed in our tests.

As we now have a way to get the model implied volatility from the local volatility (with $\Phi$), one can follow a Picard fixed-point algorithm as proposed in [38, 46] that we describe below.

Remark*.*

It is stated but not proved in [38] that the map $\Phi$ is contracting and so is $f\left(\left\{\sigma_{LV}\right\}\right)\rightarrow\left\{\sigma_{LV}\right\}*\frac{\left\{\Sigma_{Target}\right\}}{\Phi\left(\left\{\sigma_{LV}\right\}\right)}.$ Assuming this to be true, $f$ admits a unique fixed point that is the limit of the sequence of local volatility guesses $\left(\sigma_{LV}^{n}\right)_{n\in\mathbb{N}}$ defined as $\left\{\sigma_{LV}^{n+1}\right\}=f\left(\left\{\sigma_{LV}^{n}\right\}\right)$. In practice, convergence is achieved particularly fast (between 10 and 20 iterations).

The calibrated local volatility is shown in Figure E.1, where we plot it on a time scale to $T_{\max}$ for a better illustration of its shape.

We perform the calibration with 800 space steps and 100 time steps per year for the forward Dupire PDE, where we use the finite element method with quadratic basis functions. We then price quoted vanilla contracts with the backward Feynman–Kac PDE under the calibrated local volatility model. We get a maximum error in implied volatility smaller than 0.01% (i.e., for a market volatility of 20%, the calibrated volatility could be $20.00\pm 0.01$% in the worst case scenario).

Additionally, we plot the discounted marginal density of the spot extracted from the market. As mentioned before, this quantity is $\frac{\partial^{2}C_{LV}}{\partial K^{2}}$ and can be computed from the PDE solution immediately and accurately. As we will use the density in the calibration formula in Theorem 2.5, we want it to be smooth and accurate. Figure E.2 shows that quantity.

Appendix F Four-factor hybrid stochastic volatility model calibration

Consider a “purely stochastic” version of the model (2.1) – the Heston-2CIR

++

model with leverage function $\equiv 1$ – and additionally suppose that the domestic and the foreign short interest rate dynamics are independent of the dynamics of the spot FX rate. The model is governed by the following system of SDEs under the domestic risk-neutral measure $\mathbb{Q}^{d}$:

[TABLE]

where $W$ and $W^{V}$ are correlated Brownian motions with correlation coefficient $\rho$. Note that the quanto correction term in the drift of the foreign short rate vanishes due to the postulated independence assumption between the spot FX rate and foreign short rate dynamics.

Define the vector of parameters $\beta_{2}=(v_{0},\kappa,\theta,\xi,\rho)$. The next step in our calibration is to find the values of these $5$ model parameters for which European call option prices best match the market call prices retrieved from volatility quotes for different strikes and maturities. For a EURUSD transaction, the market standard is to choose USD as the domestic currency and EUR as the foreign currency. The forward FX rate for a payment date $T$ is defined as

[TABLE]

where $P^{d}(0,T)$ and $P^{f}(0,T)$ are the domestic and foreign discount factors at time [math] for a maturity $T$, respectively.

Under the postulated simple correlation structure of the Brownian drivers and when the short rates are driven by the CIR process, i.e., when $h^{d,f}=0$, Ahlip and Rutkowski [3] derive an efficient closed-form formula for the European call option price. Hence, we denote by $C_{\text{A}}(K,T)$ the fair value under the Heston–2CIR model of a European call option with strike $K$ and maturity $T$ computed with the aforementioned formula, and by $C_{\text{H}}(K,T)$ the fair value of the same option but under the Heston–2CIR

++

model. For $i\in\{d,f\}$, we define for brevity

[TABLE]

where the shift functions $h^{d,f}$ were calibrated in Appendix D. Then we can extend the pricing formula of Ahlip and Rutkowski [3] as follows.

[TABLE]

Therefore, $C_{\text{H}}(K,T)=H_{f}^{-1}C_{\text{A}}(\tilde{K},T)$, where $\tilde{K}=H_{f}H_{d}^{-1}K$. We now calibrate the Heston–2CIR

++

model by minimising the sum of the squared differences between model and market call prices:

[TABLE]

where $\sigma_{i,j}$ is the quoted volatility corresponding to a strike $K_{j}$ and a maturity $T_{i}$, for $j=1,\ldots,m$ and $i=1,\ldots,n$. There are many ways to choose the objective function (error measure) in (F.3). For instance, we may consider either call prices or Black–Scholes implied volatilities and minimise the sum of either absolute or relative (squared) differences between model and market values, using either uniform or non-uniform weights. We choose this particular error measure, which assigns more weight to more expensive options (in-the-money, long-term) and less weight to cheaper options (out-of-the-money, short-term), for two reasons. First, the Heston model, and hence the Heston–2CIR

++

model by extension, cannot reproduce the smiles or skews typically observed for short maturities that well and a more careful calibration to these smiles would result in a larger overall model error due to the inherent poor fit of the model to the short-term. Second, market data becomes scarce as the maturity increases, and hence we already assigned more weight to the short- and mid-term sections of the volatility surface; for instance, we have more maturities up to $1$ year than between $1$ and $5$ years.

As before, we employ a nonlinear least-squares solver (the trust-region-reflective algorithm, see [12]) for the calibration and a global optimisation method (a genetic algorithm) for verification purposes. Due to the non-linearity and non-convexity of the problem, the calibrated model parameters may end up in a local rather than a global minimum of the objective function. Hence, a good initial parameter guess may significantly improve the quality of the calibration. Practitioners usually use variance swap prices to calibrate $v_{0}$, $\kappa$ and $\theta$. In our case, we found the squared ATM $3$-week and $5$-year volatilities to provide good initial guesses for $v_{0}$ and $\theta$, respectively.