SLV with Stochastic Correlation & Jumps

• SLV with stochastic correlation and correlated jumps is an advanced FX modeling framework that captures time-varying skew and jump-induced discontinuities.
• The model employs coupled stochastic differential equations driven by correlated Brownian motions and Lévy jump processes to reflect dynamic market behavior.
• A fully implicit PIDE solver combined with a multi-stage calibration procedure ensures efficient and accurate pricing of vanilla and exotic FX options.

The stochastic local volatility (SLV) model with stochastic spot/volatility correlation and correlated jumps is an advanced mathematical framework for modeling the dynamics of foreign exchange (FX) option markets. The model is designed to capture stochastic skew—the empirical phenomenon that the implied volatility skew observed in FX options markets exhibits random, time-varying behavior—not reproducible by constant-correlation stochastic volatility or jump models. The approach considers spot (FX rate), instantaneous variance, and their correlation as state variables, each driven by a mixture of correlated Brownian motions and Lévy jump processes, enabling comprehensive modeling of skew and its stochastic dynamics, including jump-induced discontinuities and term-structure effects (Itkin, 2017).

1. Stochastic Local Volatility Model Formulation

The model consists of three coupled stochastic differential equations (SDEs) in risk-neutral measure:

• $S_t$: the FX spot;
• $v_t$: the instantaneous variance;
• $R_t = \operatorname{artanh}\,\rho_t$ with $\rho_t = \tanh R_t \in [–1,1]$, the spot/variance instantaneous correlation parameterized via the artanh transformation to ensure boundedness.

The SDEs are: $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$ where:

• $W_{s,t}$, $W_{v,t}$, $W_{r,t}$ are correlated Brownian motions;
• $L_{i,t}$ are pure-jump Lévy processes with both idiosyncratic and common components;
• $c, a \in [0,2)$ are power-law exponents (e.g., for Heston, $v_t$0, $v_t$1);
• $v_t$2 is a local volatility surface.

This formulation includes stochastic local volatility, mean-reverting variance and correlation, and both idiosyncratic and common jumps [(Itkin, 2017), Sect. 2].

2. Correlated Diffusion and Jump Structures

2.1 Diffusion Correlation

The Brownian drivers have instantaneous correlations: $v_t$3 where $v_t$4 is constant, ensuring the $v_t$5 instantaneous correlation matrix is positive semidefinite for all $v_t$6. The stochasticity in $v_t$7 introduces dynamic skew effects.

2.2 Correlated Lévy Jumps

Each jump process $v_t$8 is decomposed as: $v_t$9 with $R_t = \operatorname{artanh}\,\rho_t$0 independent idiosyncratic Lévy processes and $R_t = \operatorname{artanh}\,\rho_t$1 a common Lévy process, coupled by coefficients $R_t = \operatorname{artanh}\,\rho_t$2. Their pairwise jump-correlation is

$R_t = \operatorname{artanh}\,\rho_t$3

The total instantaneous correlation for each pair includes both diffusion and jump contributions, as in

$R_t = \operatorname{artanh}\,\rho_t$4

with $R_t = \operatorname{artanh}\,\rho_t$5, $R_t = \operatorname{artanh}\,\rho_t$6, and $R_t = \operatorname{artanh}\,\rho_t$7. There is no correlation between Brownian and jump parts [(Itkin, 2017), Sect. 2.2].

3. PIDE Framework for Option Pricing and Calibration

3.1 Backward PIDE (Pricing)

Option valuation is based on a three-dimensional Partial Integro-Differential Equation (PIDE) for the undiscounted option value $R_t = \operatorname{artanh}\,\rho_t$8: $R_t = \operatorname{artanh}\,\rho_t$9 where:

• $\rho_t = \tanh R_t \in [–1,1]$0 comprises convection, diffusion, and mixed-derivative terms corresponding to the SDE drift and diffusive covariance structure,
• $\rho_t = \tanh R_t \in [–1,1]$1 is an integral operator capturing jumps via the joint Lévy measure over $\rho_t = \tanh R_t \in [–1,1]$2,
• $\rho_t = \tanh R_t \in [–1,1]$3 is the discount rate.

3.2 Forward PIDE (Calibration)

For calibration, the forward density $\rho_t = \tanh R_t \in [–1,1]$4 satisfies the adjoint PIDE: $\rho_t = \tanh R_t \in [–1,1]$5 This equation allows for simultaneous calibration across strikes and maturities. The adjoint operators transpose drift, diffusion, and integral components accordingly [(Itkin, 2017), Sect. 4].

4. Fully Implicit Finite-Difference Solver

The numerical solution employs a fully implicit Strang-splitting Alternating Direction Implicit (ADI) scheme enhanced by Picard-splitting for mixed derivatives:

• Strang Split: Each time-step ($\rho_t = \tanh R_t \in [–1,1]$6) consists of half diffusion, full jumps, then half diffusion steps:
  1. $\rho_t = \tanh R_t \in [–1,1]$7
  2. $\rho_t = \tanh R_t \in [–1,1]$8
  3. $\rho_t = \tanh R_t \in [–1,1]$9

• 1D Sub-Steps: Each univariate convection-diffusion operator is solved implicitly using backward-Euler or Crank-Nicolson formalisms with $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$0 for unconditional stability and M-matrix structure. These sub-steps are linear in each spatial direction.

• Mixed Derivatives: Rather than an explicit seven-point finite-difference, the mixed second derivative $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$1 is handled by sequential 1D implicit solves, ensuring nonnegativity, second-order spatial accuracy, and stability (damping parameter $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$2, as shown in Prop. 5.2).
• Jump Operators: Split among components and executed using matrix exponentials or ADI approximations. Under a double-exponential law for jumps (e.g., Kou, 2004), the approach remains computationally efficient in $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$3 time per direction.
• Accuracy and Stability: The scheme achieves unconditional $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$4 stability, mass conservation (Preserves $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$5), and second-order accuracy in time and space by design [(Itkin, 2017), Sect. 5].

The table summarizes operator characteristics:

Operator Type	Numerical Method	Complexity/Features
1D Convection/Diffusion	Implicit (BV-Euler/CrN)	Banded, $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$6 per direction
Mixed Derivative	Picard Split	Unconditionally stable, $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$7
Jump (Lévy)	Matrix Exponential/ADI	$$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$8, negative-definite, mass-pres.

5. Multi-Stage Calibration Procedure

The model admits a four-stage sequential calibration protocol:

1. Local Volatility Surface: Calibrate $$\begin{aligned} dS_t &= (r_d - r_f) S_t\, dt + \Sigma_s(S_t, t)\, S_t^c\, \sqrt{v_t}\, dW_{s,t} + S_t\, dL_{s,t} \ dv_t &= \kappa_v(t)[\theta_v(t) - v_t]\, dt + \xi_v\, v_t^a\, dW_{v,t} + v_t\, dL_{v,t} \ dR_t &= \kappa_r(t)[\theta_r(t) - R_t]\, dt + \xi_r\, dW_{r,t} + dL_{r,t} \end{aligned}$$9 by fitting to vanilla European options with variance and correlation dynamics disabled ($W_{s,t}$0, $W_{s,t}$1, no jumps), equivalent to a Dupire local vol step.
2. SV Diffusion: With volatility-of-volatility enabled ($W_{s,t}$2), fit $W_{s,t}$3 to at-the-money options and variance swap data.
3. Correlation Diffusion: Adjust correlation parameters $W_{s,t}$4 to fit risk-reversal term structures—exotics sensitive to skew—using the forward PIDE solver for efficiency.
4. Jumps: Fit idiosyncratic jump parameters $W_{s,t}$5 for $W_{s,t}$6 to short-dated smile features; calibrate common-jump loadings $W_{s,t}$7 and variance $W_{s,t}$8 to match stochastic skew across maturities.

Optimization minimizes a weighted objective function over option prices or implied volatilities, with gradients approximated via finite differences or adjoints and solved via Levenberg-Marquardt or hybrid optimizers. The forward PIDE enables efficient all-strikes-at-once calibration [(Itkin, 2017), Sect. 5].

6. Empirical Performance and Stochastic Skew

Numerical experiments conducted on a non-uniform $W_{s,t}$9 space grid and $W_{v,t}$0 show the following qualitative and quantitative model behaviors:

• In the Heston-type limit ($W_{v,t}$1, $W_{v,t}$2), model output matches FFT benchmarks for ATM calls within 0.02% error.
• The time-varying risk-reversal skew as a function of maturity shows significant variability with stochastic correlation, including up to ±250 bps relative difference compared to constant-correlation models.
• Barrier and double-no-touch (DNT) exotic volatilities exhibit shifts of 5–50 bps depending on the stochasticity of $W_{v,t}$3.
• Introducing correlated jumps, even with low common-jump intensity ($W_{v,t}$4), produces up to ±10% shifts in barrier risk-reversal skew.

These results confirm that incorporating both stochastic correlation and correlated jumps across spot, variance, and correlation dimensions enables the reproduction of observed stochastic skew in FX options—manifested as sign changes, term-structure oscillations, and amplitude variations of the implied-vol slope at fixed delta unattainable with constant-$W_{v,t}$5 stochastic volatility plus jump models. The fully implicit, mass-conserving, Strang-ADI forward PIDE solver enables practical calibration to market data within feasible computational times ($W_{v,t}$6 minute for a full grid run on standard hardware) [(Itkin, 2017), Sect. 6].