Modeling stochastic skew of FX options using SLV models with stochastic spot/vol correlation and correlated jumps

TL;DR

This paper introduces a novel SLV model with stochastic spot/variance correlation and Levy-driven drivers to capture the stochastic behavior of FX options' implied volatility skew, supported by a stable numerical scheme.

Contribution

It develops a new SLV model with stochastic correlation driven by Levy processes and proposes a stable finite-difference scheme for calibration to market data.

Findings

Model effectively captures stochastic skew in FX options

Numerical scheme is unconditionally stable and second-order accurate

Simulations demonstrate the model's capacity to fit market prices

Abstract

It is known that the implied volatility skew of FX options demonstrates a stochastic behavior which is called stochastic skew. In this paper we create stochastic skew by assuming the spot/instantaneous variance correlation to be stochastic. Accordingly, we consider a class of SLV models with stochastic correlation where all drivers - the spot, instantaneous variance and their correlation are modeled by Levy processes. We assume all diffusion components to be fully correlated as well as all jump components. A new fully implicit splitting finite-difference scheme is proposed for solving forward PIDE which is used when calibrating the model to market prices of the FX options with different strikes and maturities. The scheme is unconditionally stable, of second order of approximation in time and space, and achieves a linear complexity in each spatial direction. The results of simulation obtained by using this model demonstrate capacity of the presented approach in modeling stochastic skew.

Tables

Table 1. Table 1: Parameters of the test for pricing an European Call option.

$T$	$r_d$	$r_f$	$L$	$H$	${\kappa }_v$	${\xi }_v$	${\theta }_v$	${\kappa }_r$	${\xi }_r$	${\theta }_r$	${\rho }_{vr}$	$S_0$	$v_0$	${\rho }_{xy,0}$
0.5	0.02	0.01	50	84.5	2	0.3	0.1	0.3	5	-0.2	0.4	65	0.5	-0.7

Equations

$$d{\rho }_t=\alpha (\overline{\rho }-{\rho }_t)+\beta \sqrt{(1-{\rho }_t)({\rho }_t-{\rho }_{min})}d{W}_t+(1-{\rho }_t)d{N}_t\cdot 1_{I_t<I_0}.$$

$$\lambda ={\lambda }_0{\left(\frac{1-I_t}{1-I_c(t)}\right)}^{\kappa },$$

$$I_c(t)=F(t)e^{-\frac{1}{2}{\sigma }^{2}t-n\sigma \sqrt{t}},$$

$$R_t=\mathrm{arctanh}({\rho }_t).$$

$$d{S}_t$$

$$d{v}_t$$

$$d{R}_t$$

$$S_{t}\Big{|}_{t=0}=S_{0},\quad v_{t}\Big{|}_{t=0}=v_{0},\quad R_{t}\Big{|}_{t=0}=R_{0}=\operatorname{arctanh}(\rho_{0}).$$

$$Af(x)={\int }_{\mathbb{R}}\left(f(x+y)-f(x)-y{1}_{|y|<1}\right)\mu (dy),$$

$${\int }_{|y|>1}{e}^y\mu (dy)<\infty .$$

$$⟨d{W}_{s,t}d{W}_{r,t}⟩$$

$$⟨d{W}_{v,t}d{W}_{r,t}⟩$$

$$⟨d{W}_{s,t}d{W}_{v,t}⟩$$

$$\mathbb{E}[W_{s,t}{W}_{r,t}]=\mathbb{E}\left[{\int }_0^t{\tilde{\rho }}_{s}ds\right]\ne {\tilde{\rho }}_{t}t$$

$$\mathcal{R}=\left(\begin{array}{ccc}1& {\tilde{\rho }}_t& {\tilde{\rho }}_{sr}\\ {\tilde{\rho }}_t& 1& {\tilde{\rho }}_{vr}\\ {\tilde{\rho }}_{sr}& {\tilde{\rho }}_{vr}& 1\end{array}\right)$$

$$det(\mathcal{R})=1-{\tilde{\rho }}_t^2-{\tilde{\rho }}_{sr}^2-{\tilde{\rho }}_{vr}^2+2{\tilde{\rho }}_t{\tilde{\rho }}_{vr}{\tilde{\rho }}_{sr}\ge 0,$$

$$-B-{\tilde{\rho }}_t{\tilde{\rho }}_{sr}\le {\tilde{\rho }}_{vr}\le B-{\tilde{\rho }}_t{\tilde{\rho }}_{sr},B=\sqrt{(1-{\tilde{\rho }}_t^2)(1-{\tilde{\rho }}_{sr}^2)}.$$

$$det(\mathcal{R})=(1-{\tilde{\rho }}_t^2)(1-{\tilde{\rho }}_{vr}^2)\ge 0,\forall ({\tilde{\rho }}_t\in [-1,1])\cup ({\tilde{\rho }}_{vr}\in [-1,1]).$$

$$X_t=(X_{1,t},...,X_{n,t}{)}^{\top }=(Y_{1,t}+b_1{Z}_t,...,Y_{n,t}+b_n{Z}_t{)}^{\top }$$

$${\varphi }_{\mathbf{X}}(\mathbf{u};t)={\varphi }_Z\left(\sum _{i=1}^n{b}_i{u}_i;t\right)\prod _{i=1}^n{\varphi }_{Y_j}(u_j;t),\mathbf{u}\in {\mathbb{R}}^n.$$

$${\rho }_{ij}=\frac{b_j{b}_i\text{mbox}Var(Z_1)}{\sqrt{\mathrm{Var}(X_{j,1})}\sqrt{\mathrm{Var}(X_{i,1})}}.$$

$${\rho }_{sv}$$

$${\rho }_{vr}$$

$${\rho }_{sr}$$

$$\mathrm{Var}(L_{i,1})$$

$$V_{\tau }=[\mathcal{D}+\mathcal{J}-r]V,$$

$$\mathcal{D}$$

$$F_1$$

$$F_2$$

$$F_3$$

$$F_0$$

$$+\left[\tanh (R){\tilde{\rho }}_{vr}\sigma (S_t)\sigma (R_t)+b_s{b}_r\mathrm{Var}(Z_1)\right]\frac{{\partial }^2}{\partial S\partial R}$$

$$+\left[{\tilde{\rho }}_{vr}\sigma (v_t)\sigma (R_t)+b_v{b}_r\mathrm{Var}(Z_1)\right]\frac{{\partial }^2}{\partial R\partial v}\equiv F_{sv}+F_{sr}+F_{vr}.$$

$$\sigma (S_t)={\Sigma }_s(S,t)S^c\sqrt{v},\sigma (v_t)={\xi }_v{v}^{\alpha },\sigma (R_t)={\xi }_r.$$

$$\mathcal{J}V$$

$$\displaystyle-\sum_{\chi\in[s,v]}(e^{y_{\chi}}-1)\frac{\partial V(x_{s},x_{v},x_{r},\tau)}{\partial x_{\chi}}-y_{r}\frac{\partial V(x_{s},x_{v},x_{r},\tau)}{\partial x_{r}}\Bigg{]}\mu(dy_{s}dy_{v}dy_{r}),$$

$$V(S,v,R,T)=P(S),$$

Full text

††journal:

1 Introduction

It is known that in the FX world skew of the implied volatility of the FX options demonstrates some kind of stochastic behavior. For instance, Carr and Wu, (2007) report that for two currency pairs (which are the U.S. dollar prices of the Japanese yen and the British pound) by analyzing the market implied volatilities of option series they find the slope of the smile to greatly vary over time. Namely, the sign of the slope switches several times in the considered sample. Therefore, although the risk-neutral distribution of the currency return exhibits persistent fat-tail behavior, the risk-neutral skewness of the distribution experiences strong time-variation which can be very positive or very negative on any given date.

To capture this effect when doing option pricing, new models are required because, as mentioned in Carr and Wu, (2007), the jump-diffusion stochastic volatility models like that in Bates, (1996), are able to capture the average shape of implied volatility smiles and the time-variation in their levels but cannot generate strong time variation in the risk-neutral skewness of currency returns. We can mention at least two approaches that are used to attack this problem.

The first one is proposed in Carr and Wu, (2004) where skew is modeled by separate the up jumps and the down jumps in spot. These jumps are modeled by two Lévy processes. The stochastic volatility and skewness are introduced by a separate time change in each Lévy component. The stochastic variation in the relative proportion of up and down jumps generates stochastic variation in the risk-neutral skewness of currency returns. The model is analytically tractable for pricing European vanilla options. However, as per Higgins, (2014) it often requires jumps in spot that are much larger than those actually experienced. Therefore, it can misprice barrier options if the probability of spot jumping through the barrier is significant Also pricing of exotic options requires numerical methods, see e.g., Itkin and Carr, (2011).

Another approach is to randomize those model parameters that govern the risk-neutral skewness. In jump-diffusion models they are the mean jump size and the correlation coefficient between the currency return and the stochastic variance process. The first parameter governs the risk-neutral skewness at short maturities, and the second one - at long maturities. Therefore, in a set of papers the correlation between the spot and instantaneous variance is considered to be stochastic and driven by a separate stochastic process, see van Emmerich, (2006); Da Fonseca et al., (2007); Ma, (2009); Ahdida and Alfonsi, (2013); Higgins, (2014); Zetocha, 2015a and references therein. Moreover, in Zetocha, 2015b where a pricing problem is considered for basket equity options, the standard correlation model, which is based on the Jacobi process, is extended by adding jumps in the correlation process. This is done based on the analysis of the market data for basket options which demonstrates that the jump to high level of correlation occurs after several down moves that drag the basket below some comfort zone of normality. The jumps are modeled in a simplified way by adding to the Jacobi process a term proportional to $dN$ with $N$ being a Poisson process whose intensity depends on the basket performance and time.

2 Modeling stochastic correlation

We mention two basic approaches to model the stochastic correlation $\rho_{t}$. Consider first only the diffusion processes. Since by definition $\rho_{t}$ is bounded in $[-1,1]$, the first approach would be to model it as a bounded diffusion. One of the popular choices is the Jacobi process used, e.g., in van Emmerich, (2006); Zetocha, 2015a . However, this requires the process to be mean-reverting, while it is not obvious why the correlation should be mean-reverted.111Perhaps, the only reason would be if one wants an expectation of the stochastic correlation to be stationary at long time horizons. However, we don’t know any paper that would justified this based on the available market data. Other bounded diffusions can also be used.

This approach is tractable for European vanilla options, but not for the exotic options. Because of that, in Higgins, (2014) an analytical approximation for barrier and one-touch options is developed which is based on semi-static vega replication. It is shown that this approximation works well in markets where the risk neutral drift is modest.

The second approach, van Emmerich, (2006), uses a transformation $\rho_{t}=g(X_{t})$ where $X_{t}$ is an arbitrary diffusion process (perhaps driftless), and function $g(X_{t})$ transforms the support of the $X_{t}$ outcomes to $[-1,1]$. Popular transformations are hyperbolic tangent, Teng et al., (2016), normal CDF, Carr, (2012), normalized inverse tangent, van Emmerich, (2006), etc. In principle, any continuous map $\mathbb{R}\to[-1,1]$ can be used for this purpose, so the real choice becomes a matter of having some additional properties, e.g., tractability, to be taken into account. Also, despite this approach is a bit less intuitive, it opens the door to various sophisticated models of the stochastic correlations.

One of such possible extensions would be a natural idea to consider the underlying process $X_{t}$ to be a jump-diffusion or a Lévy process. To the best of our knowledge, there is the only attempt to model correlation by introducing jumps which is presented in Zetocha, 2015a . However, this model is very simplistic and empirical. It is introduced based on the analysis of some market data on baskets of stocks, with an observation that the jump to high level of correlation might happen after several down moves that put the basket below some level of normality. Therefore, when modeling these baskets using the stochastic correlation model, Zetocha, 2015a defines $\rho_{t}$ to follow the extended Jacobi process

[TABLE]

Hear $\alpha,\beta,\rho_{min}$ are some constant model parameters, $W_{t}$ is the standard Brownian motion, $N_{t}$ is the Poisson jump process with the intensity $\lambda$ given by the following expression

[TABLE]

where $\lambda_{0},\kappa$ are constants, $I_{c}(t)$ is the critical level defined as $n$ standard deviations of the basket performance $I_{t}$

[TABLE]

Here $\sigma$ is the at-the-money volatility of the basket at time $t$, and $F(t)$ is the basket forward at time $t$. So the jumps can occur if the basket performance $I_{t}$ hits the level $I_{t}=I_{0}$ at $t>0$. Otherwise, if $I_{t}-I_{0}$ is positive the jump doesn’t occur. The jump size $(1-\rho_{t})$ in Eq.(1) is chosen in such a way, that if the jump occurs, $\rho_{t}$ always jumps from $\rho_{t^{-}}$ to 1 (full correlation).

Obviously, this model is empirical with a lot of parameters introduced, which financial meaning is not entirely transparent. Therefore, calibration of such a model could be a hard issue. Also the assumption that the jump always results in the full correlation to be reached immediately after the jump occurs, seems to not be justified at all.

Therefore, in this paper we utilize a different approach. We use a general Lévy models framework to introduce jumps, see e.g., Cont and Tankov, (2004). Since a general Lévy process is not bounded, we use a certain transformation as this is described in above, to bound it to $[-1,1]$. Once this is done, we consider three stochastic processes: the spot, instantaneous variance and their correlation each driven by some Lévy process. We assume that the Brownian components of each process are correlated, and the jump component are correlated similar to how this is done in Itkin and Lipton, (2015); Itkin, 2016b . However, the Brownian motions and the jump components remain uncorrelated. This model has a clear financial meaning, and can be calibrated in few steps. For instance, idiosyncratic components of jumps can be calibrated separately to the corresponding market; the local volatility functions can be calibrated to a set of FX European vanilla options; then stochastic volatility and correlation parameters (the diffusion part) of the model can be calibrated to FX exotic options which demonstrate a smile. Finally, the common jumps parameters of the model can be calibrated to the implied volatility skew of the exotic options. This procedure runs in a loop until it converges.

The main contribution of the paper consists in several results. First, a SLV model for pricing exotic FX options which demonstrate the existence of stochastic skew is proposed which assumes the spot/InV correlation to be stochastic. Moreover, the model also includes correlated Lévy jumps in each of the stochastic drivers. Second, a new fully implicit finite-difference splitting scheme is proposed to solve a forward Kolmogorov equation which is of second order in both temporal and each spatial dimension, is unconditionally stable and preserves both positivity and norm of the solution while complexity-wise is linear in the number of grid nodes in each spatial dimension. We also formulate and prove a theorem that the proposed discretization preserves the norm of the solution.

The rest of the paper is organized as follows. In Section 3 we formulate the model which considers stochastic correlation of Brownian drivers and deterministic correlation of Lévy jump processes. Section 4 discusses the corresponding backward PIDE and boundary and initial conditions. In Section 5 same consideration is given to the Forward Kolmogorov equation. In Section 6 we introduce a new fully-implicit finite-difference scheme for this forward PIDE and analyze its properties. Results of some numerical experiments obtained by using this scheme are provided in Section 7. An Appendix provides a detailed derivation of the forward FD scheme.

3 Model

We consider an LSV model with stochastic spot/InV correlation and jumps by introducing stochastic dynamics for variables $S_{t},v_{t},R_{t}$. Here $S_{t}$ is the spot price, $v_{t}$ is the instantaneous variance, and

[TABLE]

In other words, the stochastic correlation $\rho_{t}$ is represented as $\rho_{t}=\tanh(R_{t})$. Therefore, since $R_{t}\in\mathbb{R}$, this map guarantees that $\rho_{t}\in[-1,1]$.

Similar to Itkin, 2016b , we consider a model where the stochastic SDE for each variable $S_{t},v_{t},R_{t}$ includes both diffusion and jumps components, as follows:

[TABLE]

subject to the following initial conditions

[TABLE]

Here $r_{d},r_{f}$ are the domestic and foreign interest rates, $t$ is the time, $\Sigma_{s}$ is the local volatility function, $W_{s,t},W_{v,t},W_{r,t}$ are the correlated Brownian motions, $\kappa_{v},\theta_{v},\xi_{v}$ are the mean-reversion rate, mean-reversion level and volatility of volatility (vol-of-vol) for the instantaneous variance $v_{t},\kappa_{r},\theta_{r},\xi_{r}$ are the corresponding parameters for $R_{t}$, $0\leq a<2,\ 0\leq c<2$ are some constants (power constants) which are introduced to add an additional flexibility to the model as compared with the popular Heston ($\alpha=0.5$), lognormal ($\alpha=1$) and 3/2 ($\alpha=1.5$) models. As mentioned in Itkin, 2016b , one has to be careful if she wants to determine these parameters by calibration. This is because having both the vol-of-vol and the power constant in the same diffusion term brings an ambiguity into the calibration procedure. However, it can be resolved if some additional financial instruments are used in calibration, e.g., exotic option prices are combined with the variance swaps prices, see Itkin, (2013).

In the last equation of Eq.(5) the drift term is introduced to be mean-reverting. However, this might not be necessary222Or, as it was already mentioned, there is no evidence that the stochastic correlation should be mean-reverting.. Therefore, when doing so we rely on calibration of the model to the market data. In other words, if the data exhibits mean-reversion, we expect the mean-reversion level $\theta_{r}(t)$, found by the calibration, to differ from zero. On contrary, if $\theta_{r}(t)=0$, there is no mean-reversion, and the sign of the drift term depends on the sign of $\kappa_{r}(t)$.

Processes $L_{s,t},L_{v,t},L_{r,t}$ are pure discontinuous jump processes with generator $A$

[TABLE]

with $\mu(dy)$ be a Lévy measure, and

[TABLE]

At this stage, the jump measure $\mu(dx)$ is left unspecified, so all types of jumps including those with finite and infinite variation, and finite and infinite activity could be considered here.

The last line of Eq.(5) is a combination of the Hull-White model, Hull and White, (1990), for the diffusion component with a (arithmetic) Levy process for the jump component. As such, this model allows both positive and negative values of $R_{t}$.

3.1 Correlations

Below for a better transparency we consider three cases.

No jumps

Consider first the case when there is no jumps in Eq.(5). Then

[TABLE]

where we use the symbol $\widetilde{}$ to mark parameters and variables for a pure diffusion case. The last line in Eq.(7) is explained in detail in van Emmerich, (2006). It means that

[TABLE]

as it would be if $\widetilde{\rho}_{t}$ is constant.

We assume that $\widetilde{\rho}_{vr}$ is constant, but not $\widetilde{\rho}_{sr}$. Indeed, the correlation matrix

[TABLE]

must be positive-semidefinite. This provides restrictions on, at least, one of the correlation coefficients $\widetilde{\rho}_{sr},\widetilde{\rho}_{vr}$. For instance, let us assume that $\widetilde{\rho}_{sr}\in[-1,1]$ and then, since also $\widetilde{\rho}_{t}\in[-1,1]$, the semi-definiteness is defined by the condition

[TABLE]

which is solved by

[TABLE]

These allowed values of $\widetilde{\rho}_{vr}$ are depicted in Fig. 1.

So the domain of definition of $\widetilde{\rho}_{vr}$ depends on $\widetilde{\rho}_{t}$, and since $\widetilde{\rho}_{t}$ is stochastic, it also changes stochastically. Therefore, $\widetilde{\rho}_{vr}$ cannot be a constant.

To resolve this issue, following the idea of Higgins, (2014), further we assume that $\widetilde{\rho}_{vr}=const,\ \widetilde{\rho}_{sr}=\widetilde{\rho}_{t}\widetilde{\rho}_{vr}$. Thus, the spot/correlation correlation is also stochastic, but driven by the $\widetilde{\rho}_{t}$ process multiplied by a constant. This provides

[TABLE]

Pure jump processes

It is known that models with jumps in the spot price better predict the observed market data, especially at short maturities. It is also known, that jumps in $v_{t}$ might also be needed. For instance, in Sepp, 2011b ; Sepp, 2011a exponential and discrete jumps are investigated in both $S_{t}$ and $v_{t}$. The author concludes that infrequent negative jumps in $v_{t}$ are necessary to fit the market data on equity options. Definitely, jumps in both spot and the instantaneous variance could be correlated. Then it would be natural to extend this point and introduce the correlation between the jumps in $\rho_{t}$ and that in $S_{t}$, as well as between jumps in $\rho_{t}$ and that in $v_{t}$. For instance, Clements and Liao, (2013) considered the links between co-jumps within a group of large stocks, the volatility of, and correlation between their returns. They found that, despite the occurrence of common, or co-jumps between the stocks is unrelated to the level of volatility or correlation, both volatility and correlation are lower subsequent to a co-jump. Therefore, it does make sense to consider jumps in the spot, instantaneous variance and their correlation to be fully correlated.

In this paper we introduce these correlations as this is done in Itkin, 2016b ; Itkin and Lipton, (2015) by following the approach of Ballotta and Bonfiglioli, (2014). They construct the jump process as a linear combination of two independent Lévy processes representing the systematic factor and the idiosyncratic shock, respectively (see also Cont and Tankov, (2004)). It has an intuitive economic interpretation and retains nice tractability, as the multivariate characteristic function in this model is available in closed form based on the following proposition of Ballotta and Bonfiglioli, (2014):

Proposition 3.1

Let $Z_{t},\ Y_{j,t},\ j=1,...,n$ be independent Lévy processes on a probability space $(Q,F,P)$, with characteristic functions $\phi_{Z}(u;t)$ and $\phi_{Y_{j}}(u;t)$, for $j=1,...,n$ respectively. Then, for $b_{j}\in\mathbb{R},\ j=1,...,n$

[TABLE]

is a Lévy process on $\mathbb{R}^{n}$. The resulting characteristic function is

[TABLE]

By construction every factor $X_{i,t},i=1,...,n$ includes a common factor $Z_{t}$. Therefore, all components $X_{i,t},i=1,...,n$ could jump together, and loading factors $b_{i}$ determine the magnitude (intensity) of the jump in $X_{i,t}$ due to the jump in $Z_{t}$. Thus, all components of the multivariate Lévy process $\mathbf{X}_{t}$ are dependent, and their pairwise correlation is given by (again see Ballotta and Bonfiglioli, (2014) and references therein)

[TABLE]

Such a construction has multiple advantages, namely:

As $\mbox{sign}(\rho_{i,j})=\mbox{sign}(b_{i}b_{j})$, both positive and negative correlations can be accommodated 2. 2.

In the limiting case $b_{i}\to 0$ or $b_{j}\to 0$ or $\mathrm{Var}(Z_{1})=0$ the margins become independent, and $\rho_{i,j}=0$. The other limit $b_{i}\to\infty$ or $b_{j}\to\infty$ represents a full positive correlation case, so $\rho_{i,j}=1$. Accordingly, $b_{i}\to\infty,\ b_{3-i}\to\infty,\ i=1,2$ represents a full negative correlation case as in this limit $\rho_{i,j}=-1$.

If we want the jump correlations to be constant, then this setting is sufficient. However, since we want the correlation $\rho_{sv}\equiv\rho_{t}$ to be stochastic, this approach has to be modified. In this paper we don’t elaborate on it, and use the former one, so assuming that the correlation between jumps is constant.

General case

When the underlying Lévy processes have both the diffusion and jumps components, the correlation matrix can be build as a combination of two constructions presented in above. We assume that all Brownian motions $W_{i,t},\ i\in[S,v,r]$ are independent of $L_{i,t}$. In this setting the total instantaneous correlations read

[TABLE]

where $\sigma_{s},\sigma_{v},\sigma_{r}$ are volatilities of the diffusion processes of $S_{t},v_{t},R_{t}$, and $Y_{i,t}$ are the idiosyncratic jump processes. With such a definition all correlations $\rho_{vr},\rho_{sr},\rho_{t}$ could take any value within the interval $[-1,1]$ with no restrictions. As follows from Eq.(11) $\rho_{vr}$ is constant while $\rho_{t},\rho_{sr}$ are stochastic. If the second term in the numerator of $\rho_{sv}$ in Eq.(11) is much bigger that the first one, this correlation also is constant and determined by the correlation of jumps. If, on the contrary, the second term is small as compared with the first one, $\rho_{sv}$ is stochastic and defined by the correlation of the Brownian motions. A similar argument is valid for $\rho_{sr}$ as well.

4 Pricing PIDE

To price contingent claims, e.g., vanilla or exotic options written on the underlying spot price, a multidimensional (backward) PIDE could be derived by using a standard technique as in Cont and Tankov, (2004). It describes the evolution of the option price $V$ under risk-neutral measure and reads

[TABLE]

where $\tau=T-t$ is the backward time, $T$ is the time to the contract expiration, $r$ is the risk free interest rate, $\mathcal{D}$ is the three-dimensional linear convection-diffusion operator of the form

[TABLE]

[TABLE]

According to Eq.(11)

[TABLE]

In Eq.(12) $\mathcal{J}$ is the jump operator

[TABLE]

where $\mu(dy_{s}dy_{v}dy_{r})$ is the three-dimensional Lévy measure, and $x_{s}=\log S/S_{0}$, $x_{v}=\log v/v_{0},\ x_{r}=R$.

This PIDE has to be solved subject to the boundary and terminal conditions. Consider first a pure diffusion case (no jumps). For the FX European vanilla options the terminal condition reads

[TABLE]

where $P(S)$ is the vanilla option payoff as defined by the corresponding contract (Call or Put). For the double barrier options considered in this paper, the payoff reads

[TABLE]

where $\tau_{H}$ is the first time the spot price hits the upper barrier $H$, and $\tau_{L}$ is the first time the spot price hits the lower barrier $L$. There also exist contracts where some rebate is paid when the spot price hits either the upper or lower barriers. So Eq.(16) can be easily extended to this case.

The boundary conditions also depend on the type of the options under consideration. In this paper we consider European vanilla and double barrier options. The boundary conditions for those differ only by the conditions at the boundaries of the spot price domain. Namely, we set

[TABLE]

for the double barrier options, and

[TABLE]

for the vanilla Call options.

For the instantaneous variance $v_{t}$ the form of the boundary conditions depends on the exponent $a$ in Eq.(5). For instance, when $a=1/2$, i.e. the diffusion part of $v$ is the CIR process, it is known that no boundary condition is required at $v=0$ if the Feller condition $2\kappa_{v}\theta_{v}>\xi_{v}^{2}$ holds. So in this case the equation Eq.(12) with $v=0$ is used as the boundary condition. More general, if $a,\kappa_{v},\xi_{v},\theta_{v}$ are chosen in such a way that the process $v$ never reaches the origin, no boundary condition is required. Otherwise, either the reflection or absorbing boundary conditions can be set at $v=0$. At $v\to\infty$ the boundary condition is $\partial V/\partial v\to 0$. Alternatively, it can be written $V(S,v,R,t)\|_{v\to\infty}=S$. In more detail, see, e.g., discussions in Itkin and Carr, (2011); Haentjens and In’t Hout, (2012); de Graaf, (2012).

In the presence of jumps these conditions should be extended as follows. Suppose we want to use finite-difference method to solve the above PIDE and construct a jump grid, which is a superset of the finite-difference grid used to solve the diffusion equation (i.e. when $\mathcal{J}=0$, see Itkin, 2016a ). Then these boundary conditions should be set on this jump grid as well as at the boundaries of the diffusion domain.

Boundary conditions for $R_{t}$ should be set at $R\to\pm\infty$ which corresponds to $\rho_{t}\to 1$ and $\rho_{t}\to-1$. In the pure diffusion case (no jumps) at $\widetilde{\rho}_{t}=-1$ from Eq.(9) we obtain $\widetilde{\rho}_{sr}=\widetilde{\rho}_{vr}=0$. Thus, in words, this is the case when the spot price is perfectly anti-correlated with the instantaneous variance, while both are independent of their stochastic correlation, which is 1. As follows from Eq.(5), the spot price dynamics then is defined by a local volatility model, Dupire, (1994); Derman and Kani, (1994), with the local volatility function $\Sigma_{s}(S_{t},t)S_{t}^{c}\sqrt{v_{t}}$. Therefore, no boundary condition is necessary at $R\to-\infty$, since the PDE in Eq.(12) will itself serves as the boundary condition at this end. At $\widetilde{\rho}_{t}=1$, we have $\widetilde{\rho}_{sr}=\rho_{vr}$. Again, no boundary condition is necessary at $R\to\infty$, and the PDE in Eq.(12) itself serve as the boundary condition with those values of correlations substituted into it.

When the jumps are taken into account, perfect correlation $\rho_{l}=1$ is achieved when the loading factors $b_{s},b_{v}$ both tend to infinity. Hence, contributions to Eq.(11) from both the diffusion correlation and from the idiosyncratic parts are very small. It then follows from Eq.(11) that

[TABLE]

In addition, all terms in the definition of $F_{0}$ in Eq.(13) tend to infinity when $b_{s}\to\infty,b_{v}\to\infty$. Therefore, we must set

[TABLE]

Hence, similar to the pure diffusion case, we can assume that $V=V(S,\tau)$ at $R\to\infty$. Therefore, Eq.(12) in this limit again transforms to that for the local volatility model, and this equation is just the boundary condition at $R\to\infty$.

By a similar argument, $V=V(S,\tau)$ at $R\to-\infty$. Therefore, substituting this into the PIDE, we obtain the equation which is the boundary condition at this end.

5 Forward PIDE

It is well-known that the backward PDE/PIDE is helpful for pricing options, e.g., with the same strike and maturity but various initial spot prices. However, when calibrating some model to the options market data, we need to simultaneously price multiple options with different strikes and maturities but the same initial spot price. For doing so, a forward Kolmogorov equation can be used. Forward PIDEs can be derived using techniques proposed in Cont and Bentata, (2012); Andersen and Andreasen, (2000); Carr, (2006). Alternatively, we can exploit the approach which is discussed in Itkin, (2015) and for the financial literature is originated by Lipton, (2001, 2002), see a literature survey in Itkin, (2015). Briefly, the idea is as follow.

Instead of a continuous forward equation we consider a discretized Kolmogorov equation from the very beginning. In other words, this is equivalent to the discrete Markov chain defined at space states corresponding to some finite-difference grid $\mathbf{G}$. Consider first only the stochastic processes with no jumps. It is known (see, e.g., Goodman, (2004)) that the forward Kolmogorov equation for the density of the underlying process $X(t)$ can be written in the form

[TABLE]

where $p({\bf s},t)$ is a discrete density function (i.e., a vector of probabilities that $X(t)$ is in the corresponding state ${\bf s}$ at time $t$), $T$ is the expiration time, and the generator $\mathcal{A}$ has a discrete representation on $\mathbf{G}$ as a matrix of transition probabilities between the states. For a given interval $[t,t+\Delta t]$ where the generator $\mathcal{A}$ does not depend on time, the solution of the forward equation is

[TABLE]

Therefore, to compute the option price $V=e^{-rT}E_{Q}[P(X(T))]$, where $P(X(t))$ is the payoff at time $t$ and $E_{Q}$ is an expectation under the risk-neutral measure $Q$, we start with a given $p(0)$, evaluate $p(T)$ by solving Eq.(19),333If $\mathcal{A}=\mathcal{A}(t)$, we solve in multiple steps in time $\Delta t$ by using, for instance, a piecewise linear approximation of $\mathcal{A}$ and the solution in Eq.(20). and finally evaluate $e^{-rT}\int_{0}^{\infty}p(s,T)P(s,T)ds$.

On the other hand, one can start with a backward Kolmogorov equation for the option price

[TABLE]

which by change of variables $\tau=T-t$ could be transformed to

[TABLE]

It is also well-known that $\mathcal{B}=\mathcal{A}^{\top}$; see, e.g., Goodman, (2004).

Based on that, constructing the transposed discrete operator $\mathcal{A}^{\top}$ on a grid $\mathbf{G}$ is fairly straightforward. However, extending this idea to modern finite-difference schemes of higher-order approximation (see, e.g., In’t Hout and Welfert, (2007) and references therein) is a greater obstacle. This obstacle is especially challenging because these schemes by nature do multiple fractional steps to accomplish the final solution, and because an explicit form of the generator $\mathcal{A}$, that is obtained by applying all the steps, is not obvious. Finally, including jumps in consideration makes this problem even harder.

This problem is solved in Itkin, (2015) where a splitting finite-difference scheme for the 2D forward equation is constructed which is fully consistent in the above-mentioned sense with the backward counterpart. For the latter, two popular scheme are chosen, namely: Hundsdorfer and Verwer (HV) and a modified Craig-Sneyd scheme, In’t Hout and Welfert, (2007). Also a consistent forward scheme is proposed in Itkin, (2015) for solving the corresponding PIDEs.

6 Finite-difference scheme for a forward PIDE

To solve the forward PIDE, we again use the approach of Itkin, (2015). It consists of few steps

6.1 Splitting of diffusion and jumps

We use the Strang splitting scheme, Strang, (1968) which provides the second-order of approximation in a temporal step. For the backward equation it reads 444In practical implementations, when solving the backward equation, the term $-rV$ is usually included into the operator $\mathcal{D}$ by splitting it into three equal parts $-rV/3,$ and adding each part to $F_{1},F_{2}$ and $F_{3}$, correspondingly.

[TABLE]

which could be re-written as a set of fractional steps:

[TABLE]

Thus, instead of an unsteady PIDE, we obtain one PIDE with no drift and diffusion (the second equation in Eq.(24)) and two unsteady PDEs (the first and third ones in Eq.(24)).

To produce a consistent discrete forward equation, we use the transposed evolutionary operator, which results in the scheme

[TABLE]

assuming that $\Delta\tau=\Delta t$. The fractional steps representation of this scheme is

[TABLE]

6.2 Splitting scheme for the pure diffusion equation

Here we consider the solution of the equation

[TABLE]

which is required at the first and third steps of Eq.(26). For doing that in Itkin, (2015) two forward finite-difference schemes are proposed for the 2D case. The first scheme, which is consistent with the HV scheme is

[TABLE]

where the subindex $n$ marks the $n$-th step in time, $M_{i}\equiv I-\theta\Delta tF^{n}_{i}$, $F_{i}^{n}=F_{i}(t_{n})$, and $\Delta\widetilde{Y}_{i}=\widetilde{Y}_{i}-Y_{i-1},\ i=1,2$, and $\theta$ is the parameter of the HV scheme.

The forward analog for another popular backward finite-difference scheme—a modified Craig-Sneyd (MCS) scheme, see In’t Hout and Foulon, (2010), can also be found in Itkin, (2015).

Despite both these schemes are unconditionally stable, there exists a delicate issue regarding the positivity of the solution. As mentioned in Itkin, (2015), for steps 1, 2, and 4 in Eq.(28) this is guaranteed if both $M_{1}^{\top}$ and $M_{2}^{\top}$ are M-matrices; see Berman and Plemmons, (1994). To achieve this, an appropriate (upward) approximation of the drift (convection) term has to be chosen, which is often discussed in the literature; see In’t Hout and Foulon, (2010) and references therein.

For steps 3 and 5, a possible issue is that the central difference approximation of the second order for the mixed derivatives doesn’t preserves positivity of the solution. Various attempts are made in the literature to replace this approximation with a reduced scheme which uses a 7 point stencil, again see the discussion in Itkin, (2015). Lack of positivity usually gives rise to instability unless a very small spatial step is used which is impractical. This is less important for the backward scheme because the explicit step is followed by some number of implicit steps which significantly damp possible instabilities (so the resulting solution is non-negative). But for the forward scheme the explicit step is the last one, and so this could be critical.

To resolve this problem, in Itkin, 2016b it is proposed to replace the explicit step for the mixed derivatives with the implicit one. It was proved there that the scheme proposed for calculating of the mixed derivative terms guarantees the positivity of the solution, while the whole scheme becomes fully implicit. Below we describe this scheme for the backward equation while for the forward one it is exactly same if one changes the notation from $V$ to $p$, and from $\tau$ to $t$, and also replaces all matrices with their transposes.

6.2.1 Fully implicit scheme for the mixed derivatives

As applied to a 3D case, the idea is to represent the explicit step of the HV scheme Eq.(56) as

[TABLE]

where $F=\sum_{j}F_{j},\ j\in[0,3]\cap\mathbb{Z}$. Now observe, that the rhs of Eq.(29) is a Padé approximation (0,1) of the equation

[TABLE]

The solution of this equation can be formally written as

[TABLE]

Alternatively, a Padé approximation (1,0) can also be applied to all exponentials in Eq.(31) providing same order of approximation in $\Delta\tau$ but making all steps implicit. Namely, this results into the following splitting scheme of the solution of Eq.(30):

[TABLE]

Steps 2-4 in Eq.(32)) are the one-dimensional implicit steps, so we already know how to solve these equations.

By construction of the HV scheme, the first step in Eq.(32) has to be solved with the accuracy $O(\Delta\tau)$, again see discussion in Itkin, 2016b . With this accuracy, the previous equation can be factorized as

[TABLE]

or using splitting

[TABLE]

The order of the splitting steps usually doesn’t matter.

For illustration, let us consider a mixed derivative term $F_{Sv}$ in Eq.(33). We start with replacing the explicit step used in the HV and MCS schemes for the mixed derivative term

[TABLE]

by the implicit representation, which has the same order of approximation (first order in time)

[TABLE]

where $W(S)=\sqrt{\mathrm{Var}(S_{t})},W(v)=\sqrt{\mathrm{Var}(v_{t})}$.

We re-write this equation in the form

[TABLE]

where $P,Q,$ are some positive numbers which have to be chosen based on some conditions, e.g., to provide diagonal dominance of the matrices in the parentheses in the lhs of Eq.(36).

As shown in Itkin, 2016b , Eq.(36) can be solved using fixed-point Picard iterations. One can start with setting $V(\tau+\Delta\tau)=V(\tau)$ in the rhs of Eq.(36), then solve sequentially two systems of equations

[TABLE]

Here $V^{k}$ is the value of $V(\tau+\Delta\tau)$ at the $k$-th iteration.

The following main theorem provides a necessary discretization of the above equation.

Proposition 6.1** (Proposition 2.3 from Itkin, 2016b )**

Let us assume $\rho_{s,v}\geq 0$, and approximate the lhs of Eq.(6.2.1) using the following finite-difference scheme:

[TABLE]

Then this scheme is unconditionally stable in time step $\Delta\tau$, approximates Eq.(6.2.1) with $O(\max(h^{2}_{S},h^{2}_{v}))$ and preserves positivity of the vector $V(x,\tau)$ if $Q=\beta\sqrt{\Delta\tau}/h_{v},\ P=\beta\sqrt{\Delta\tau}/h_{S}$, where $h_{v},h_{S}$ are the grid space steps correspondingly in $v$ and $S$ directions, and the coefficient $\beta$ must be chosen to obey the condition:

[TABLE]

The scheme Eq.(38) has a linear complexity in each direction.

For the proof, see Appendix B in Itkin, 2016b . Here we define a one-sided second order approximations to the first derivatives: backward approximation $A^{B}_{2,x}:\ A^{B}_{2,x}V(x)=[3V(x)-4V(x-h)+V(x-2h)]/(2h)$, and forward approximation $A^{F}_{2,x}:\ A^{F}_{2,x}V(x)=[-3V(x)+4V(x+h)-V(x+2h)]/(2h)$. Also $I_{x}$ denotes a unit matrix. All these definitions assume that we work on a uniform grid with step $h$, however this could be easily generalized for the non-uniform grid as well, see, e.g., In’t Hout and Foulon, (2010).

A practical choice of the coefficient $\beta$ is discussed in Itkin, 2016b . It is also shown there that in case $\rho_{s,v}<0$ the same Theorem holds, if in the terms containing $\rho_{s,v}$ we replace $A^{F}_{2,S}$ with $A^{B}_{2,S}$ and vice versa.

6.2.2 Presence of jumps

When jumps are taken into account, and hence the loading factors $b_{i}\neq 0,\ i\in[s,v,r]$, according to Eq.(13) all mixed derivatives operators will include an extra term. For instance, Eq.(34) will now read

[TABLE]

where $F_{Sv}(\tau)$ represents the diffusion part of the covariance. Thus, this case is reduced to the previous one with no jumps, since by assumption, $b_{i}$ and $\mathrm{Var}(Z_{1})$ are deterministic.

6.2.3 Fully implicit scheme for the forward equation

Using this idea, we derive a forward scheme for the 3D diffusion equation which is an exact adjoint of the backward HV scheme. This is done similar to the derivation in Itkin, 2016a for the 2D case. However, in the present approach we combine it with the fully implicit scheme for the mixed derivative term. The derivation of the implicit-explicit scheme is given in Appendix. The final result reads

[TABLE]

where $\Delta Y=\widetilde{Y}-Y$

To apply Proposition 6.1, we rewrite the fourth line of Eq.(40) in the form

[TABLE]

Since we are constructing a scheme which provides the second order of approximation in $\Delta\tau$, the result of the multiplication $\overline{Y}_{i}=L^{n}_{i}Y_{i}$ is equivalent up to $O((\Delta t)^{2})$ to the solution of the equation

[TABLE]

see discussions in Itkin, 2016b about Padé approximations and properties of M-matrices. The Eq.(42) is a pure implicit equation, and its solution is a non-negative vector provided that $M_{i}^{\top}$ is an M-matrix, and vector $Y_{i}$ is also non-negative.

Thus, Eq.(41) can be re-written as

[TABLE]

The equation for $Z_{0}$ takes a form of Eq.(29). Therefore, Proposition 6.1 can be applied to compute the vector $\widetilde{Y}_{0}$ in a fully implicit manner which preserves the positivity of the solution, and is of second order of approximation in the spatial steps.

In a similar way we can proceed with the last line of Eq.(40) which can be represented as

[TABLE]

Introduce $\widehat{Y}_{i}$ as the solution of the equation

[TABLE]

Observe, that since the solution of Eq.(44) has to be obtained with the accuracy up to $O((\Delta t)^{2})$, we can replace operators $L^{n-1}_{i}$ in the second sum of the last line of Eq.(44) by $L^{n}_{i}$, because $\Delta tF^{n-1}_{i}=\Delta tF^{n}_{i}+O((\Delta t)^{2})$. For the same reason $Z_{1}=Z_{0}+O((\Delta t)^{2})$. Then

[TABLE]

The last term in this equation, again has the form of Eq.(29) and can also be computed using the fully implicit scheme described in Proposition 6.1, if one replaces $A_{2}^{B},A_{2}^{T}$ by their transposes. This finalizes the construction of the fully implicit splitting scheme for the 3D advection-diffusion equation which provides second order of approximation in temporal and spatial steps, is unconditionally stable and preserves positivity of the solution.

As shown in Itkin, 2016b , this fully implicit scheme allows elimination of first few Rannacher steps as this is usually done in the literature to provide a better stability (see survey, e.g., in Haentjens and In’t Hout, (2012)), and provides much better stability of the whole scheme which is important when solving multidimensional problems.

For the sake of convenience, below we collect all the above splitting steps just in one block

[TABLE]

Complexity-wise, this scheme requires the solution of 12 implicit systems of linear equations, and 2 explicit equations. The latter steps being treated using the fully-implicit scheme require per one explicit equation: 3 implicit steps for $i=1,2,3$, approximately 2 iterations for each of 3 mixed terms, and for each term solving 2 implicit systems, so totally 15 implicit systems. Overall, the whole solution could be obtained by solving 42 implicit systems555By re-grouping terms in Eq.(47) it actually can be reduced to 36.. Therefore, the complexity of this approach is about $42\cdot O(N_{1}N_{2}N_{3})$ where $N_{i}$ is the number of the grid nodes in the $i$-th direction, $i=1,2,3$. So this is about 2.3 times slower than the original HV scheme for the backward equation.

6.3 Preserving the norm of the forward solution

Since the discrete solution of the forward Kolmogorov equation is the discrete probability density, it has to obey two important conditions. The first one is the positivity of the solution and was considered in the previous section. The other condition requires that the three-dimensional integral of the density should be equal to 1. In the discrete case this can be replaced with the condition that

[TABLE]

Then the natural question would be: do there exist any conditions that an appropriate finite-difference scheme should obey in order to automatically preserve Eq.(48)? In the one-dimensional case the answer to this question is given by the following Proposition:

Proposition 6.2

Consider a forward scheme Eq.(47). A typical splitting step of this scheme consists in solving an equation of the type

[TABLE]

where $Y_{0}$ is a given vector, $Y$ is the unknown vector to be determined, and $M=I-\theta\Delta tF$ is the lhs matrix. The vector $Y_{0}$ has all positive elements, and its $L_{1}$ norm is 1, i.e.

[TABLE]

In other words, $Y_{0}$ belongs to the probability vectors $\mathbb{P}$, $Y_{0}\in\mathbb{P}$, or to the discrete distributions. Then $Y$ is also a discrete distribution, i.e. $Y\in\mathbb{P}$, if discretization of $F$ is such that the matrix $B=\|F\|$ is a valid generator matrix, $B\in\mathbb{G}$ (see Itkin, (2017) and references therein), and $M$ is an EM- (or an M-)666Every M-matrix is also an EM-matrix. matrix.

• Proof

Since $B\in\mathbb{G}$, it has the following property of the generator matrix, Itkin, (2017)

[TABLE]

If $B$ satisfies Eq.(49), then obviously

[TABLE]

On the other hand, one can represent $Y$ as $Y=(M^{\top})^{-1}Y_{0}$. Therefore, to have $Y$ with the same $L_{1}$ norm as that of $Y_{0}$, the matrix $(M^{\top})^{-1}$ should have its columns to be in $\mathbb{P}$, $(M^{\top})^{-1}_{:,j}\in\mathbb{P}\ \forall j\in[1,N]$.

Observe, that since $M$ is an EM-matrix, $M^{\top}$ is also an EM-matrix, and $(M^{\top})^{-1}$ is a positive matrix, again see Itkin, (2017). Now the result that $(M^{\top})^{-1}_{:,j}\in\mathbb{P}\ \forall j\in[1,N]$ follows from the fact that if $Y_{0}\in\mathbb{P}$, $(M^{\top})^{-1}Y_{0}$ is a convex combination of the columns of $(M^{\top})^{-1}$ with coefficients given by the entries of $Y_{0}$. Each column of $(M^{\top})^{-1}$ must be in $\mathbb{P}$, and $\mathbb{P}$ is a convex set.

Further we need the following Lemma:

Lemma 6.3

Suppose we are given an invertible matrix $M^{\top}$ with real entries with the sum of every column is 1. Then, for the matrix $(M^{\top})^{-1}$ the sum of the elements in each column is 1.

• Proof

The statement of this Lemma is obviously equivalent to the following statement. Given an invertible matrix $M$ with real entries with the sum of every row is 1, for the matrix $M^{-1}$ the sum of the elements in each row is also 1. To prove this, denote $M^{-1}=D$. Then we have

[TABLE]

and thus

[TABLE]

Therefore, since by construction $M^{\top}$ is such a matrix that its every column sums to 1, based on the above Lemma it follows that $(M^{\top})^{-1}\in\mathbb{P}$. This proves the statement of the Proposition. $\blacksquare$

In 3D case the splitting scheme in Eq.(47) consists of multiple 1D steps. Therefore, intuitively, it is clear that if all discretizations of matrices $F_{i},\ i=0,1,2,3$ are chosen according to the Proposition 6.2, the total scheme preserves the $L_{1}$ norm of the solution vector. A rigorous prove of this statement will be given elsewhere. This is also confirmed by our numerical experiments.

6.4 Splitting scheme for the pure jump equation

For the second (jump) step in Eq.(24) a splitting scheme is proposed in Itkin and Lipton, (2015); Itkin, 2016b . If, e.g., jumps in $S$ are described by an exponential Lévy process, the jump integral admits representation in the form of a pseudo-differential operator

[TABLE]

which is introduced in Itkin and Carr, (2012), and where $x=\log S$. In the definition of the operator $\mathcal{J}$ (which is actually an infinitesimal generator of the jump process), the integral can be formally computed under some mild assumptions about existence and convergence if one treats the term $\partial/\partial x$ as a constant. It is important to emphasize that

[TABLE]

where $\psi(u)$ is the characteristic exponent of the jump process, and MGF(u) is the moment generation function corresponding to this characteristic exponent. This directly follows from the Lévy-Khinchine theorem. This point of view apparently has been pioneered in Jacob, (1996). See also the detailed state-of-the-art surveys in Jacob and Schilling, (2001); Böttcher et al., (2014). Then, using Proposition 3.1, the operator $e^{\Delta\tau\mathcal{J}}$ in Eq.(24) can be represented as, Itkin and Lipton, (2015); Itkin, 2016b

[TABLE]

Thus, this requires a sequential solution of 7 equations at every time step.

As shown in Itkin, (2015), using this method for the forward equation doesn’t bring any problem, because

[TABLE]

and there is no issue with computing $\mathcal{J}^{\top}$. Indeed, if $\mathcal{J}$ is the negative of an M-matrix, the transposed matrix preserves the same property. Also if $\mathcal{J}$ has all negative eigenvalues, the same is true for $\mathcal{J}^{\top}$. Then the unconditional stability of the scheme and its property of preserving the positivity of the solution follow from Proposition 4.1 in Itkin, (2013). Also, it is easy to show that the jump equation at the central step of Eq.(53)

[TABLE]

can be solved with respect to the unknown vector $Z$ by using the same ADI method as in Itkin and Lipton, (2015); Itkin, 2016b , by simply replacing all matrices $\|\triangledown_{\eta}\|$ with $\|\triangledown_{\eta}\|^{\top},\ \eta\in[s,v,r]$.

6.5 Boundary and initial conditions

As mentioned in Itkin, (2015), the boundary conditions for the forward scheme should be consistent with those for the backward scheme. However, these two sets of conditions are not exactly same because for the forward equation the dependent variable is the density, while for the backward equation the dependent variable is the undiscounted option price.

In Section 4 we have already discussed the boundary and initial (terminal) conditions for the backward equation. For the forward equation the obvious initial condition reads

[TABLE]

where $\delta(x)$ is the Dirac delta-function.

Setting the boundary conditions is a more delicate issue. In the $S$ domain, we require the density $p(S,v,R,t)$ to vanish at the domain boundaries, i.e. at $S=0$ and $S\to\infty$. The boundary conditions in the $v$-domain are discussed, e.g., in Lucic, (2008). For the Heston model, if the Feller condition $2\kappa_{v}\theta_{v}>\xi_{v}$ holds, the boundary condition at $v=0$ is $\partial p(S,v,R,t)/\partial v=0$. Otherwise, $p(t,S,0,R)=0$ should be used as the boundary condition. Therefore, in this paper we use this approach, i.e. if the point $v=0$ is inaccessible based on given values of the parameters $a,\kappa_{v},\xi_{v},\theta_{v}$, we use the boundary condition $\partial p(S,v,R,t)/\partial v=0$ at $v=0$. At $v\to\infty$ we set $p(S,v,R,t)\to 0$.

Also, since by construction the $R$ domain doesn’t contain any reflection or absorbing boundaries, the analysis of the option price behavior at $R\to\pm\infty$ given in Section 4 can be applied here as well. Based on that, we conclude that similar to the $S$ domain, we can require the density $p(S,v,R,t)$ to vanish at $R\to\pm\infty$.

7 Numerical experiments

In this section we describe some results of our numerical experiments. Here, for the sake of brevity we provide just a single example where jumps are taken into account, all the other examples deal with a pure diffusion case. We will present more detailed results for the model with jumps elsewhere.

Fo all experiments we solve the forward Kolmogorov equation Eq.(19) using the proposed finite-difference fully implicit splitting scheme, and compute the density $p(S,v,R.t)$. Parameters of the model used in these tests are given in Tab. 1, where $L,H$ denote the lower and upper barriers, correspondingly. The local volatility function $\sigma(S,t)$ in these tests is always set to 1. Power parameters are taken to be $a=0.5,c=1$. Thus, this setting is analogous to the Heston model, but with the stochastic correlation.

The finite-difference non-uniform grid is constructed similar to how this was done in Itkin, 2016b and includes 101 nodes in $S$, and 81 nodes in $v$ and $R$ directions, so the whole grid is $101\times 81\times 81$. The fixed temporal step is 0.01. No Rannacher or any other smoothing scheme is used at the first temporal steps. Parameters $\beta$ for computing each mixed derivatives term by using the fully implicit scheme are taken to be 10. Parameter $\theta$ of the HV scheme according to Itkin, (2015) is 0.3. We use Matlab and run our code at PC with Dual Quad Core Intel(R) Core(TM)i7-4790 CPU 3.60 Ghz.

Since the proposed scheme is constructed by using the transposed operators (matrices), and a fully implicit scheme for the mixed derivatives, its convergence analysis coincides with that provided in Itkin, 2016b . Therefore, in our numerical experiments we are mostly dedicated to the financial interpretation of the results obtained.

7.1 European options

In this tests we applied our model to pricing European vanilla options. As a benchmark the Heston model is chosen which is mimicked by setting $\kappa_{r}=\xi_{r}=\rho_{yz}=0$ in the general setting described in above.

The ATM Call option value in this test at $T=0.5$ computed using the forward equation is 10.7533 while the benchmark value computed using FFT is 10.7564. Thus, for this test the relative accuracy is about 0.02%. An elapsed time for computing one step in time is about 4.8 sec.

Fig. LABEL:eurSkew presents the RR(10) skew of the European options computed using the proposed model with parameters as in Table 1 (the SS1 model) and the Heston model. It can be seen that the Heston model also generates stochastic skew. This is the known fact,777I appreciate our discussion with Peter Carr on the subject. because the skew is proportional to the spot/InV covariance. Therefore, even if the correlation coefficient is constant, the skew changes with time if the spot variance, or the vol-of-vol, or both change with time (so their product changes with time).

It is also seen that making the correlation stochastic doesn’t contribute much to the magnitude of the skew. In Fig. LABEL:eurSkewDif the difference in skew between the Heston and SS1 models is presented together with the third case which displays the difference in skew between the Heston model and the stochastic skew model (SS2) where in contrast to Table 1 we use $\kappa_{r}=3,\theta_{r}=0.,\xi_{r}=5,\rho_{yz}=0.4$. In other words, in the SS2 case there is no mean-reversion in the dynamics of the correlation coeffcient $\rho_{t}$. The relative difference in skew as compared with the Heston model could reach +100% or -250%, however, this is because the absolute value of the skew is small.

Fig. LABEL:Eur3DXY,LABEL:Eur3DXZ,LABEL:Eur3DYZ display the computed density of the European vanilla options at maturity. Since the whole picture is 4D, the results are presented in various planes. In other words, the third dimension is represented as a sequence of graphs corresponding to some discrete values of the third independent variable. It can be seen, that the density in various planes could be bimodal, or even quadro-modal. Such type of the solutions for the density function has been already observed, e.g., in Ashyraliyev et al., (2008); Kumar and Narayanan, (2006).

7.2 Down-and-out barrier options

It is known, that when the lower barrier $L$ is fixed, for the DO barrier Call option there always exists a strike corresponding to the option $\delta$ = 0.1 (10 Delta). However, for the Put this stops to be true when $T$ increases. So starting from some maturity $T$, the equation $\partial V(S,K,v,\rho,T,r_{d},r_{f})/\partial S=0.1$ has no solution with respect to $K$. Therefore, in practice to trade the skew the barrier options are often set with the barriers being a function of time, Kainth and Saravanamuttu, (2016).

In the next test we set a linear time-dependent lower barrier to be $L(t)=L-20t$, where $L$ is given in Table 1. The corresponding results for the skew of the down-and-out barrier options are presented in Fig. LABEL:doSkew,LABEL:doSkewDif. Here, we again use a non-uniform grid with 100 x 80 x 80 nodes. Projections of the grid in planes $S-v$ and $S-R$ are presented in Fig. LABEL:gridDO.

It can be seen that assuming the correlation coefficient $\rho_{sv}$ to be stochastic bring some changes in the skew value as compared with the pure stochastic volatility Heston model. At the short time scale these changes are within $\pm$50 bps. Therefore, if the market demonstrates larger variations of the skew, adding jumps to the correlation could make the model more flexible.

7.3 Double no-touch options

In this test we compute prices of double no-touch (DNT) options with fixed the lower and upper barriers which values are given in Tab. 1. Again, we use a non-uniform grid 100 x 80 x 80 which projections in planes $S-v$ and $S-R$ are presented in Fig. LABEL:gridDNT.

In Fig. LABEL:dntIV the implied volatility of the DNT options is presented as a function of time to maturity $T$ computed for the Heston and stochastic correlation models. At the given time scale the difference between two models is about 5-10 bps and also changes the sign with time.

7.4 Down-and-out barrier options (with jumps)

This test uses the setting of Section 7.2, but now includes jumps into consideration. We represent all jumps - common and idiosyncratic - by using the Kou double exponential model, Kou and Wang, (2004), similar to how this was done in Itkin and Lipton, (2015).

The Lévy density of the double-exponential jumps is, Cont and Tankov, (2004):

[TABLE]

where $\varphi$ is the jumps intensity, $\theta_{1}>1$, $\theta_{2}>0$, $1>p>0$; the first condition was imposed to ensure that the underlying asset price has a finite expectation.

Accordingly, for this process $\mathrm{Var}(L_{i,1})$ in Eq.(11) in the explicit form reads

[TABLE]

Parameters of the jumps model for this particular test are presented in Table LABEL:TabJumpParam

For the diffusion part we again use a non-uniform grid 101 x 81 x 81. The jump non-uniform grid in each direction is a superset of the diffusion grid built using a geometric progression, see Itkin, 2016b . So the whole jump grid includes $117\times 105\times 103$. The fixed temporal step is 0.01. The typical elapsed time obtained at the same computer is: the first diffusion step - 4.7 secs, the jumps steps in $S,v,R$ - $S$ - 0.15, 0.25, 0.4 secs, the common jumps step - 5.5 secs to converge to the relative accuracy $10^{-2}$, the next jump steps in $S,v,R$ - 0.40, 0.48, 0.58 secs, and the last diffusion step - 4.7 secs. Thus, it takes about 3.5 times more to compute the solution at one step in time as compared with the case with no jumps. This coincides with the theoretical estimation of Itkin and Lipton, (2015).

The ADI scheme for computing the common and idiosyncratic jumps for the Kou model is described in Itkin and Lipton, (2015). The value of the ADI parameter $s$ is experimentally chosen to be $s=10000$ which is sufficient to provides convergence of the scheme.

The result of simulation are presented in Fig. LABEL:doSkewJump where SSJ stays for the stochastic skew model with jumps. It can be seen that even jumps with relatively small intensity as in Table LABEL:TabJumpParam change the skew by about 10%. Obviously, the proposed model has the capacity to produce lower of higher skews by varying the model parameters.

8 Conclusion

The paper deals with an advanced model which is build on top of an LSV model by including the following components: i) the spot/instantaneous variance correlation is stochastic, ii) all stochastic drives - spot, instantaneous variance and their correlation are modeled by Lévy processes, iii) all diffusion components are correlated as well as all jump components, however there is no correlation between diffusion and jumps.

This model is used to price FX options and replicate stochastic behavior of the options’ skew which was observed in market data. For the purpose of fast calibration, we consider a forward PIDE and propose a new fully implicit splitting finite-difference scheme to solve it which is fully consistent with the corresponding solution of the backward PIDE. The scheme is unconditionally stable, of second order of approximation in time and space, and achieves a linear complexity in each spatial direction. It also preserves the positivity and norm of the density function which is proven by a series of Propositions. All these results are new.

The results of simulation obtained by using this model demonstrate capacity of the presented approach in modeling stochastic skew. However, it is worth mentioning that despite the proposed model is sufficiently rich, it contains many free parameters. Therefore, calibration of the model, if done at once, could be time-consuming even when using the forward scheme. However, if various market data is available, it is better to calibrate various pieces of the model separately, as this was discussed, e.g., in Ballotta and Bonfiglioli, (2014). Namely, the idiosyncratic jumps first can be calibrated to some marginal distributions using the appropriate instruments. Then the parameters of the common jumps can be calibrated to the option prices, while keeping parameters of the idiosyncratic jumps fixed. Also the LSV part can be calibrated to the vanilla and exotic option prices as this is usually done, Bergomi, (2016).

Acknowledgments

I thank Peter Carr for various useful comments and discussions and Igor Halperin for some helpful remarks. I assume full responsibility for any remaining errors.