A new class of multidimensional Wishart-based hybrid models

In this article, we present a new class of pricing models that extend the application of Wishart processes to the so-called stochastic local volatility (or hybrid) pricing paradigm. This approach combines the advantages of local and stochastic volatility models. Despite the growing interest on the topic, however, it seems that no particular attention has been paid to the use of multidimensional specifications for the stochastic volatility component. Our work tries to fill the gap: we introduce two hybrid models in which the stochastic volatility dynamics is described by means of a Wishart process. The proposed parametrizations not only preserve the desirable features of existing Wishart-based models but significantly enhance the ability of reproducing market prices of vanilla options.

of quoted options. However, in practice we only observe prices for a discrete set of instruments with different combinations of moneyness and time to maturity. We then need to rely on techniques able to "complete" the market surface consistently with the absence of arbitrage opportunities (roughly speaking this would mean that the corresponding local variance surface remains nonnegative). Some recent papers go in this direction: in Gatheral and Jacquier (2014), easy to check conditions are formulated to eliminate static arbitrages in the famous stochastic volatility inspired (SVI) parametrization (Gatheral 2004), while a finite difference method is used in Andreasen and Huge (2011) to get a complete surface of arbitrage consistent option values from sparse market prices. Despite the ability in matching the spot implied volatilities, the LV approach has severe limitations in pricing exotic options. The inherent deterministic nature of volatility does not allow to satisfactorily describe the underlying asset volatility dynamics. A well-known consequence is that LV models tend to generate forward implied volatility surfaces much flatter than the market ones (see Gatheral 2011 for an overview on the topic).
The drawbacks of LV are partially overcome by the use of SV models where an additional stochastic process is introduced to describe the dynamics of volatility. A massive research activity has been produced in this field revealing the wide diffusion of such approach both within practitioners and academics. Notwithstanding the parametrization adopted, the more realistic description of volatility in SV models compared to the LV case comes at a price: (diffusion) processes with embedded SV usually are not able to match market implied volatilities in deep far from the money regions for short dated options. The first attempt to improve the ability of SV models to match marginal distribution of stock prices has been the use of jumps (Bates 1996). Another solution could originate from the observation that LV and SV models are complementary to some extent: this is exactly the idea behind the so-called stochastic local volatility (SLV) or hybrid approach. Combining both LV and SV models could lead to an overall improvement of pricing accuracy since the resulting hybrid model would benefit from a satisfactory description of volatility dynamics (SV) and the ability to reproduce market plain vanilla option prices (LV). A fast growing interest in SLV models has recently led to a widespread research activity aimed at formulating valid approaches. We refer to Cui et al. (2018) and Homescu (2014) for a comprehensive overview of the topic. An important feature of SLV models is that it is possible to separately calibrate the LV and SV components and combine them in a nonlinear and nonparametric way by means of the so-called leverage function. The peculiar and most delicate step to perform when dealing with SLV is the calibration of leverage function itself since we need to account for the dependency between LV and SV. Several techniques have been devised to tackle the problem: in Engelmann et al. (2011) and Ren et al. (2007) algorithms are presented to numerically solve the corresponding nonlinear integral Fokker-Planck PDE. These methods are well suited for models with embedded 1-factor SV component while for multifactor ones they would suffer from the curse of dimensionality. Recently some simulation-based techniques have been proposed: Markovian projection (Henry-Labordere 2009), particle method (Cozma et al. 2019;Henry-Labordere 2012, 2013) and nonparametric numerical scheme (Van der Stoep et al. 2014). These last two methods allow to deal easily with multidimensional processes since they basically require only to be able to sample trajectories for the asset price.
A common SLV specification in the single-asset framework is the use of nonparametric local volatility combined with a Heston-like stochastic volatility process as in Tian et al. (2015) and Van der Stoep et al. (2014). However, in the light of the poor ability of 1-factor SV models in dealing with skew dynamics, as stated in La Bua and Marazzina (2019), it seems appropriate to introduce a richer structure for the SV part. To reach such a goal, we present a new hybrid model, the Wishart stochastic local volatility model (WSLV).
Moreover, the pricing of derivatives with a multi-asset risk exposure requires a satisfactory modelling of the correlation structure among the involved assets, consistently with an accurate reproduction of individual plain vanilla market evidences. This leads to the introduction of multi-asset SLV models as done in De Col and Kuppinger (2014), where a simplified multi-asset Heston framework is extended to accommodate for hybrid asset covariations. In order to incorporate a more flexible dependence structure among state variables, we also introduce the Wishart Stochastic Local Covariance model (WSLC) as hybrid generalization of the pure SV model in Da Fonseca et al. (2007). As far as we know, this new class of models represents one of the first attempts, and probably the most comprehensive one, to provide a multidimensional variance dynamics within the SLV pricing paradigm.
Besides the extension of SLV literature to the case of matrix-based stochastic volatility parametrizations, we aim at improving the pricing performance of pure SV Wishart models. As shown in La Marazzina (2019, 2021), indeed, we must introduce stringent parameters restrictions to satisfy existence and uniqueness conditions for the solution of Wishart stochastic differential equation (SDE). This is needed, for example, when we want to simulate the variance process in order to price more exotic derivatives. However, such conditions are not usually met when market calibration is performed. Consequently, the constrained parameters set obtained enforcing such conditions is not able to accurately reproduce the market implied volatility surface. In our new setting, the additional LV component acts as a compensator meant to reduce the gap with the market still preserving the well definiteness of Wishart processes. Therefore, the introduction of the LV component in this hybrid paradigm results in a significant improvement of the ability of the model in reproducing the market.
The article is organized as follows: in Sects. 2 and 3, we briefly describe the general SLV and Wishart volatility frameworks, respectively. In Sect. 4, we present the WSLV model and report a realistic implementation of the model using the market dataset of DAX option prices. In Sect. 5, we illustrate the main properties of the WSLC model and show an application on market data for the pair of indices EuroStoxx50-DAX. Finally, in Sect. 6, we state some concluding remarks.

The general SLV framework
A SLV model specifies the following risk-neutral dynamics for the T -forward price f (t) of an equity asset where V (t) is a (possibly multifactor) stochastic process that affects the dynamics of f (t) through the functional form ψ(·) and B(t) is a Brownian motion of suitable dimension that can be correlated with the source(s) of randomness in V (t). As an example if V (t) is described in terms of 1-factor CIR process and ψ(v) = √ v, then we get the popular Heston stochastic local volatility model (Tian et al. 2015; Van der Stoep et al. 2014).
In this paper, we consider, without loss of generality, risk-free rate and dividend yield equal to zero. We refer the reader to Guyon and Henry-Labordere (2011) for the extension of SLV models to the case of stochastic interest rates and/or discrete dividends.
Here, σ (t, f (t)) is the leverage function that can be interpreted as a (deterministic) local adjustment factor introduced to match the underlying asset terminal distribution observed on the market. The definition of σ (t, f (t)) is the key aspect of SLV modelling as well as the most challenging task in the implementation of any model of this kind. Before stating the fundamental result about SLV modelling, we recall an important theorem in stochastic calculus that represents the main tool to develop the formal ground of the framework.
Theorem 1 (Mimicking process, Gyöngy 1986) Let ξ(t) be a stochastic process satisfying the SDE where W (t) is a (possibly multidimensional) Brownian motion and α 0 , α 1 bounded, nonanticipative processes such that α 1 α 1 is uniformly positive semidefinite. Then there exists a SDE with nonrandom coefficients that admits a weak solution x(t) having the same onedimensional distribution as ξ(t) for any t ≥ 0. The coefficients a 0 and a 1 have the following representation Thanks to this result, we are supplied with an explicit strategy to construct stochastic processes driven by deterministic coefficients that "mimic" some other (more involved) processes with stochastic coefficients. We now state a general result, firstly shown in Ren et al. (2007), that links the leverage function to both the LV and SV components.
Proposition 1 (Leverage function, Ren et al. 2007) Given a SLV model defined by Eq. (1), the market observed implied volatility surface is matched if and only if the leverage function is of the form 1 where σ LV (t, f ) is the local volatility calibrated on the same market prices.
Proof In its basic intuition, Eq.
(2) is a straightforward application of Theorem 1: let us set such that the solution of LV process d f (1) share the same marginal distribution for f (t). The matching with the market follows directly from the definition of local volatility process as given in Dupire (1994). We refer the interested reader to the original paper (Ren et al. 2007) for a complete proof.
Proposition 1 provides an interesting interpretation for the leverage function as the ratio between the local volatility and the expected value of the stochastic volatility component. Intuitively this suggests how a SLV model works: if the (expected) level of stochastic volatility is significantly different to the local volatility, then the leverage function acts as compensator (on average). Clearly, the specification (1) includes the "pure" LV or SV models as particular cases: with σ (t, f ) = 1 we get a standard SV model while assuming no randomness in V (t) we end up with the LV dynamics.
As pointed out in Guyon and Henry-Labordere (2011), substituting Eq. (2) into Eq. (1) we get a McKean SDE where the diffusion coefficient depends not only on state variables but also on their joint probability distribution. Keeping on following (Guyon and Henry-Labordere 2011), issues about existence and uniqueness of SLV SDE (1) have been raised since we have no certainty at all that any market implied volatility surface can be matched given a set of parameter for the SV part. Numerical investigations show that the problem becomes quite relevant when large values of volatility of volatility are needed to obtain very steep forward skew in standard (1factor) SV models. The same observation is found in Sepp (2011) where jumps are added to the SLV dynamics.

Calibration of a SLV model
In order to be able to implement a given SLV model, we firstly need to calibrate its three main components: local volatility surface, stochastic volatility parameters and leverage function (2). Remarkably, it is possible to decompose the overall problem into three individual tasks to be tackled separately. In particular, the calibration of leverage function can be carried out via ad hoc algorithms after having determined the LV and SV components. In the following, we illustrate the calibration of LV surface and leverage function. The estimate of SV parameters is, indeed, highly influenced by modelling assumptions and requires a customized approach. In our case, we refer to La Marazzina (2019, 2021) for an extensive treatment of the methodologies applied to the calibration of Wishart-based parameters.

Calibration of LV surface
The LV surface can be estimated following one of the (numerous) approaches proposed in literature. Despite the lack of evidences in favour of a particular calibration scheme, we remark that the chosen approach should not generally affect the outputs of SLV modelling, provided that the method produces realistic results. We choose to consider a nonparametric specification of LV surface and calibrate it using the algorithm in Abasto et al. (2013). The scheme, which is basically a slight modification of Andreasen and Huge (2011), has been proven to give accurate results in repricing calibration instruments with both PDE and Monte Carlo methods. Let 0 = T 0 < T 1 < T N be a set of option expiry dates and consider the Dupire forward equation (Dupire 1994) where we set, for the sake of simplicity, interest rates and dividends to zero. In Andreasen and Huge (2011), the authors propose to calibrate Eq. (4) between two consecutive dates with the following one-step implicit finite difference scheme whereσ i (K ) are the volatility proxies to be calibrated by minimizing the difference between the market prices and those obtained by Eq. (5). Further, the functionσ (T , K ) is considered to be piecewise constant in time and strike. The modification proposed in Abasto et al. (2013) considers the same optimization problem but introduces a finer time grid between market expirations in order to produce more accurate results. The corresponding equation becomes with N T the number of points in the temporal mesh between any two consecutive dates. An additional advantage of the method (in comparison with the original scheme in Andreasen and Huge 2011) is that the calibrated volatility proxies represent in this case the calibration output with no need of further processing. In Fig. 1 we report the calibrated LV surfaces for DAX and EuroStoxx50 indices. In both cases, results are obtained in less than 2 seconds via Matlab code on a laptop PC with an Intel Core i7 CPU and 8 GB RAM.

Calibration of leverage function
Given a LV surface and a set of SV parameters (ideally) calibrated on the same basket of market options, the last step in the implementation of a SLV model is the calibration of Eq.
(2). The task, however, is not trivial since the computation of the conditional expectation therein would require to know the joint distribution of f (t) and V (t). This is clearly not the case due to the presence of the nonparametric LV component that affects the dynamics of f (t) itself. In 1-factor SLV models PDE methods are proposed to solve the corresponding Fokker-Planck equation. However, in multifactor SLV models these methods are totally unfeasible since we would deal with high-dimensional PDEs.
Simulation methods appear to be more appropriate and flexible enough to be applied within multifactor SLV frameworks. Among the proposed approaches, the nonparametric method presented in Van der Stoep et al. (2014) is the most intuitive and easiest to implement. This approach is a particular implementation of the particle method with a simple rectangular kernel (Guyon and Henry-Labordere 2013). Since we extensively use this technique for our numerical tests, a brief description is given. We refer the interest reader to the original paper for a detailed illustration and a complete error analysis.
Let f be the discrete time approximation to the SLV (1). A naive Euler scheme would take the form for the generic j-th simulation, j = 1, . . . , N , at a given time step m + 1 with m = 0, . . . , M T . Here, N is the total number of simulated paths and M T the number of time steps. We also set a uniform time grid = T /M T where T is the terminal date. The (multidimensional) standard normal variable B is assumed to be enriched with the correlation structure assumed by the chosen SV model (at this stage, for sake of generality, we leave it unspecified). The stochastic variance process V is simulated with an appropriate scheme as well. Substituting Eq.
(2) in Eq. (9), we get The idea in Van der Stoep et al. (2014) . The approximation of conditional expectation in Eq. (10) is then given by for a proper choice of bins Two different bins specifications are proposed in Van der Stoep et al. (2014): equidistant and equally weighted (i.e. into each subset there is approximately the same number of realizations) bins. The latter has been found to give better results in terms of approximation accuracy in particular and it is the one used in our numerical tests. This means that for each time step we choose the + 1 bins according to indicates the sequence of f j (t m ) sorted in ascending order.

The Wishart volatility framework
In this section, we briefly introduce the Wishart (multidimensional) stochastic volatility model (WSV, La Bua and Marazzina 2019 and references therein) and the Wishart affine stochastic correlation model (WASC, La Bua and Marazzina 2021 and references therein). These two models exploit a Wishart process to deal with a multifactor stochastic volatility process in a single-asset framework (WSV) and with a multi-asset stochastic volatility framework (WASC). In Sects. 4 and 5, we show how local volatility can be used to extend the WSV and the WASC model, and the great advantages due to this extension.

Wishart process
First of all, we report the definition of a Wishart process.
Brownian motions) and S + d (R) the set of real d × d positive semidefinite matrices. We define the Wishart process as the solution on S + d (R) of the following SDE: In order to embed mean-reversion and stationarity, we consider matrix M to have only eigenvalues with negative real part. Furthermore, we relate the deterministic part of the drift in Eq. (13), , to the expected long-term value of the process, denoted with ∞ , by means of the equation From Eqs. (13) and (14), we can easily see the close connection between Wishart and CIR processes. If we set d = 1 in Eq. (13), we end up with a scalar CIR process defined by the SDE with κ, θ , and η strictly positive parameters, v 0 ≥ 0 and w v (t) a scalar Brownian motion. Section 2.1 in La Bua and Marazzina (2019) deals with the existence and uniqueness of the solution of the Wishart process SDE above defined. More precisely, assuming a more restrictive parametrization for the deterministic part of the drift there exists a unique (weak) solution if and a unique (strong) solution if The real positive parameter β plays the role of Feller's condition in the univariate case. Additionally if the Condition (17) is not met the whole process is not well defined. We refer to La Bua and Marazzina (2019) for further details.

Wishart multidimensional stochastic volatility model
In order to overcome inherent limitations of 1-factor SV models in describing the term structure of volatility skew, as documented for example in Christoffersen et al. (2009), Cont and Da Fonseca (2002) and Pacati et al. (2018), a matrix generalization of Heston model was proposed in Da Fonseca et al. (2008). This is the Wishart (multidimensional) stochastic volatility model (WSV) where the dynamics of the forward price of an equity asset is given by and with Z (t) another matrix of Brownian motions independent of W (t) and R ∈ M d (R) that fulfils the condition I d − R R ∈ S + d (R). This correlation structure is required in order to preserve the analytical tractability of the model: the WSV model so defined, indeed, belongs to the class of affine models (La Bua and Marazzina 2019).
In Da Fonseca et al. (2008), the authors show that the WSV model can be expressed in a scalar form, as stated in the following proposition.
Proposition 2 (Scalar version of WSV dynamics, Proposition 3.2 in La Bua and Marazzina 2019) Let y(t) = log ( f (t)) be the asset log-price, then its dynamics in the WSV model can be written as where b(t) and w(t) are two scalar Brownian motions with stochastic correlation given by

Wishart affine stochastic correlation model
The Wishart Affine Stochastic Correlation (WASC) model was introduced in Da Fonseca et al. (2007) and further studied in La Bua and Marazzina (2021) with the purpose of reproducing well-known multi-asset stylized facts in a tractable way. The WASC model makes use of Wishart process to describe the stochastic variance covariance matrix of asset returns. The model proposes the following joint dynamics for a vector of forward asset prices: with z(t) another vector Brownian motion independent on W (t) and r ∈ [−1, 1] d such that r r ≤ 1. Here, r can be interpreted as the vector of coefficients meant to drive the linear correlation between the shocks on asset returns and shocks on variancecovariance matrix (t). The choice of the correlation structure (25) represents the major improvement with respect to the model in Gourieroux and Sufana (2004) and aims at accommodating realistic single-asset volatility skews still preserving the affinity of the model. Remarkably, the resulting WASC dynamics (24) allows for stochastic correlation among asset returns in a tractable framework where each asset is enriched with a stochastic volatility behaviour consistent with the effects observed on plain vanilla markets. We refer to La Bua and Marazzina (2021) for further details.

The Wishart stochastic local volatility model
In La Bua and Marazzina (2019), the authors provide a fast and accurate algorithm to calibrate the WSV model. However, numerical results setting d = 2 shows that the Condition (17) is not usually met calibrating the model (e.g. Table 1 in La Bua and Marazzina 2019), and imposing it as a constrained has a not negligible impact on the ability of the model to capture the market prices and volatility (see Figures 5 and 6 in La Bua and Marazzina 2019). However, if Condition (17) is not satisfied, the process is not well defined. Therefore, in the following we exploit the LV term in order to enhance the WSV model, showing that, thanks to the LV component, we are able to fit market data also imposing Condition (17).

Model dynamics
In order to increase the ability of WSV model in describing the marginal probability distribution of asset price, we introduce the Wishart stochastic local volatility (WSLV) model given by where (t) is a Wishart process described by Eq. (13). Equivalently, we can write with dynamics of V(t) and stochastic correlation given, respectively, by Eqs. (22) and (20). The latter specification leads to an immediate interpretation of Eq. (27) as a SLV model in light of the general framework defined by Eq. (1). That is we get the WSLV dynamics as soon as we set Consistently with Proposition 1, we define the WSLV leverage function as such that, by Theorem 1, the model matches the market implied probability distribution of asset prices. The resulting model induces a hybrid instantaneous variance for the log-asset y(t) = log( f (t)) given by with the remarkable property to embed a multifactor specification of the stochastic volatility component. Here, the notation [·, ·] refers as usual to the quadratic covariation of two stochastic processes. Equation (26) represents, indeed, a generalization of the WSV model meant to combine the flexibility granted by the underlying Wishart process with an accurate pricing of plain vanilla options. Even if our modelling choice can be seen questionable since we lose the analytical tractability of WSV model we have numerous arguments in favour. Firstly, as in any SLV model, the calibration of LV and SV is performed separately, therefore we can still rely of fast calibration routines for the SV parameters. Secondly Monte Carlo pricing is unavoidable when we deal with exotic options even in affine models. Finally, as shown in Sect. 4.2, we get an evident improvement in the pricing performance of European options that could not be achieved otherwise. We think indeed that moving towards a hybrid framework is the only viable way to improve the accuracy of the WSV model given the restrictions on the parameter β.
The usual fix of adding jumps presents several drawbacks in this case. The additional increase in the number of parameters comes with no guarantee that the loss of accuracy imposed by the restriction β ≥ d − 1 would be generally overcome. Analogously to the 1-factor SV case, apart from anything else, jumps pose nontrivial problems in a risk management perspective. As illustrated in Da Fonseca and Grasselli (2011) adding jumps on the stock and/or on the volatility could lead to a loss of parameters sensitivity with respect to the skew since it is into some extent transferred to the parameters that drive the discontinuities.

Numerical results
Pricing routine within the WSLV model follows the same steps as in any other SLV model: (1) calibration of LV and SV components, (2) calibration of leverage function and (3) pricing. In La Bua and Marazzina (2019), the authors performed the calibration of Wishart-like SV component over a set of market prices of call options written on DAX index. In particular, we set d = 2 and we consider here the parameters obtained enforcing the condition β ≥ 1 (rightmost column of Table 1, reporting the parameters presented in Table 1 in La Bua and Marazzina 2019) in order to deal with a well-defined variance process, given that Condition (17) is satisfied. Exploiting the simulation-based approach in Van der Stoep et al. (2014) for the calibration of leverage function, we can address steps (2) and (3) in a single run enhancing the computational efficiency of the procedure. In order to discretize Eq. (27) we can apply the GVA sampling scheme provided in La Bua and Marazzina (2019). The use of an efficient scheme, indeed, is greatly recommended in this case: since in the SLV framework we deal with the simulation of a process with an embedded local volatility component, which is known to require a fine time discretization to get precise results, we need to rely on a fast and accurate sampling procedure. To take into account the presence of the additional local adjustment term, we need to introduce some slight modifications to the GVA scheme, that for the WSLV model reads as We are finally ready to test the ability of WSLV model in improving the pricing performance of both its components. We can therefore price European and forward starting options to assess the overall effect of the proposed SLV mixing.

European options
We test the WSLV performance with respect to short and long dated European call options. Figure 2 shows the results for maturities in the range {1M, 3M, 1Y , 3Y } in order to assess the performance over the entire term structure of volatility surface. As comparison we also priced the options within the pure LV and SV models. For the latter, we consider both parameters sets reported in Table 1, with and without the constraint β ≥ 1. In Section A, we report the numerical results in terms of difference between market and models implied volatilities. Table 2 summarizes the evidences found in terms of mean absolute error with respect to market implied volatilities: The conclusion is that if we enhance the WSV model with a LV term, then the worsening induced by the constraint on the parameter β in the WSV model can be fully neutralized: the marginal asset probability distribution in WSLV substantially matches the market one. 3 We see that the performance of WSLV is comparable to the pure LV model. This is particularly evident if we focus on the shorter maturities: for options with strike prices far from the at-the-money level, the improvement over the SV parametrizations is quite substantial. The result is a direct consequence of the SLV mechanism: the compensation effect due to the leverage function fills the gap between market and SV model. For longer maturities, since the SV model already gives adequate results, the adjustment due to the local volatility term is less pronounced. However, also in this case the WSLV gives a better fit of market values than the WSV with β ≥ 1.
We would like to stress that we also report results for the WSV model without the constraint β ≥ d − 1 to show that, without this constraint, the model performs better, as also shown in La Bua and Marazzina (2019). However, if β ∈ [0, 1), even if we can use the characteristic function for calibration, the model is not well defined, and, for example, a Monte Carlo simulation is not feasible. This is the reason why in the WSLV model we only deal with the case β ≥ 1. Our results show that, thanks to the LV component, a good fit to market data is possible even in a framework where the Wishart SV component is well defined. Therefore, introducing the LV component is a simple and efficient way to improve the fit to market data in a Wishart SV (well defined) framework.
Computationally the affine nature of WSV makes the pricing of European options quite fast without having to rely on simulations. However, the possibility to exploit new efficient simulation scheme for Wishart-based models allows to get prices in acceptable time also within the WSLV one. CPU time is then comparable to the ones for the pure SV case: the additional calculation time due to the presence of the local volatility term is indeed negligible.

Forward starting options
Having confirmed the ability of WSLV model in overcoming the limitations of the pure WSV one with respect to the pricing of European claims, we can now test if the new framework is also able to produce realistic forward volatility dynamics. More precisely, we price two different sets of forward starting options that starts, respectively, at T 1 = 3 months and T 1 = 2 years with maturity in one year. We can exploit the analytical tractability of WSV model to price forward starting options via FFT in the pure SV case. As shown in Da Fonseca et al. (2008), the forward characteristic function of log-returns can be expressed in terms of the characteristic function of the Wishart process. This gives us an insight about the nature of this kind of options: they are basically pure variance derivatives. On the other hand, we need to rely on Monte Carlo simulations for the pricing with LV and WSLV models. Figure 3 illustrates the results in terms of Black-Scholes implied volatility, that is the value of σ BS such that where C F S Model (K , t 0 , T 1 , T 2 ) is the model price of a forward starting option that starts at T 1 and expires at T 2 with strike K while C BS (S 0 , K , T , σ ) is the Black-Scholes price of a plain vanilla call option. The message we get is clear: within the WSLV framework, we are still able to reproduce the proper volatility dynamics of the pure SV model. As we can see in Fig. 3, as time goes by the forward implied volatility in the In both cases, T 2 = T 1 + 1Y LV case becomes flatter showing the inadequacy of the model in describing how the volatility actually evolves. This does not happen for the SV and WSLV models where the skew persists through time. The ability of preserving a forward implied volatility shape similar to the SV case is a well-known feature of SLV models, as shown for example in Engelmann et al. (2011) and Van der Stoep et al. (2014) for the Heston stochastic local volatility model. The problem seems poorly addressed in the existing literature, see, e.g. (Guyon and Henry-Labordere 2013, Section 12.11): according to our knowledge, this is the first time that such a result is obtained for a multifactor affine specification of the SV component. Figure 2 shows that, when an European call option is considered, the WSLV and LV models provide similar implied volatility curves, the same for the SV model, with the exception of short maturities in the case the constraint β ≥ 1, necessary for a well-defined process, is considered. In the case of forward starting call options, Fig.  3 shows that both WSV and WSLV models provide similar implied volatility curves, while LV model fails in reproducing the volatility smile, in particular if large values of T 1 are considered. To get more insight on this, we deal with the following numerical experiment: instead of calibrating the SV model on market data, we consider the four sets of parameters reported in Table 3, and we calibrate only the LV component (and the leverage function) of the WSLV model to the implied volatility market surface of European options; then we obtain the forward start (with maturity one year) implied volatility surfaces represented in Fig. 4. These results show that the implied volatilities are very different for large values of T 1 , therefore the calibration of the LV component itself is not enough to represent correctly option prices (with the exception of European   Table 3. T 1 = 3M (left), T 1 = 2Y (right). In both cases, T 2 = T 1 + 1Y call options, since numerical results here not reported show a good fit also if only the LV component of the WSLV model is calibrated). To gain more insights, we perform similar experiments considering set IV in Table 3, and comparing the LV, SV and WSLV models pricing forward start call options with 2 years maturity. Figure 5 shows the strong action of the LV component (and of the leverage function) of the WSLV model to "correct" the not-calibrated SV component, resulting on a smile which is close to the LV model. However, again, if T 1 is large, the WSLV succeeds in avoiding the flattening of the LV model. The above experiments show the importance of calibrating both the SV and the LV part correctly. In our calibrations, we calibrate both components on European options, however, even if rather standard, this is not the unique approach. One can calibrate the SV part on different kinds of derivatives (if a fast and accurate pricing method is available, e.g. if the price can be computed exploiting the characteristic function) leaving to the LV component the fit to the European implied volatility.

On the calibration of the SV and LV components of the WSLV model
These results can be generalized for any SLV model, e.g. the Heston stochastic local volatility model in Van der Stoep et al. (2014). Numerical experiments here not reported show a limited gain in exploiting the WSLV model with respect to the Heston stochastic local volatility model if the goal is to deal with derivatives like forward starting options. We would like to stress that the main advantage in exploiting matrix stochastic volatility model is the stochastic correlation between assets and volatility: indeed, in a multidimensional Heston model, the correlation is stochastic but depends on factors that generates the volatility dynamics itself, while in the case of the Wishart volatility model, the correlation depends on the volatility factors but also on other factors, see, e.g. (Benabid et al. 2008). Therefore, the advantage of using a Wishart SV component is its ability to price correctly derivatives like variance swaps, as shown in Da Fonseca et al. (2015). Due to the presence of closed formulas for some variance  Table 3. T 1 = 3M (left), T 1 = 2Y (right). In both cases, T 2 = T 1 + 2Y derivatives, see, e.g. Da Fonseca et al. (2015), in the Wishart SV framework, the SV parameters could be set to improve the fitting to variance derivatives market data, if available.

The Wishart Stochastic Local Covariance model
In La Bua and Marazzina (2021), the authors showed that the reduced WASC model is able to generalize the Heston model in a multi-asset framework embedding a rich dependence structure among all the state variables involved. Unfortunately, in order to grant a realistic dynamics of Wishart generated cross-asset correlations, they found necessary to enforce the (very) restrictive condition of positive definiteness for the Wishart process, exactly as in the WSV model in La Bua and Marazzina (2019). This in turn has the unpleasant effect to produce an extremely rigid implied volatility structure. However, even in this case, the introduction of a LV term is able to fix the problem.

Model dynamics
With the objective of alleviating such a distortion, we now introduce the Wishart Stochastic Local Covariance (WSLC) model, where the risk-neutral joint behaviour of a d-dimensional vector of forward asset prices is with (t) and b(t) defined, respectively, by Eqs. (13) and (25). Here, the matrix function [f(t), (t)] combines the local and the stochastic volatility components such that the instantaneous covariance matrix among log-assets is with y(t) = log ( f 1 (t), f 2 (t), . . . , f d (t)) . A direct application of Theorem 1 guarantees that the single-asset terminal probability distributions induced by Eq. (31) are consistent with market ones provided that we have for i = 1, . . . , d and σ LV ,i is the local volatility calibrated to the ith asset option prices. In order to fulfil this requirement, we assume the following parametrization for : where σ (t, f(t)) is a vector-valued leverage function whose values are given by It is worthwhile to point out that the choice of is not unique and different specifications can accommodate alternative modelling features. If we consider, for example, the FX market, a desirable property of a pricing model would be to preserve the symmetry under inversion and triangulation of FX rates. In that case, matrix in Eq. (34) is not the optimal choice since it would lead to mis-specified cross-rates dynamics. We refer to De Col and Kuppinger (2014) for the treatment of the topic in a multi-asset hybrid extension of (a simplified) Heston model. Nonetheless, since our main concern here is to focus on equity markets, parametrization (34) represents an adequate setting. In particular, considering for simplicity d = 2, we have that the WSLC asset instantaneous covariance matrix reads as where we use the shorthand notation σ i = σ i t, e y i (t) . The asset covariance is then hybrid (thus the name of the model), since it now depends on y 1 (t), y 2 (t) and 12 (t). Furthermore, quite interestingly, we end up with the same correlation structure induced by the pure stochastic volatility model. From Eq. (36), indeed, it results that This is due to the peculiar form of matrix . An inherent advantage of this result is that we are granted with cross-asset correlations that lie in the interval [−1, 1] as soon as the Wishart process is well defined. For the sake of generality, though, we remark that the overall effects can be more subtle and would require a thorough analysis of the induced dependence structure that we leave to further research. In particular, it would be interesting to study the differences in terms of pricing of multi-asset securities with respect to hybrid specifications of asset correlation. The framework proposed naturally applies to the most general parametrization of Wishart process. However, we consider a reduced form by setting matrix M diagonal. In this setting, the WSLC model is able to match the prices of univariate path-dependent options generated with a hybrid Heston model whose parameters are set as illustrated in Section 3.1 in La Bua and Marazzina (2021). Indeed, as already pointed out for the WASC model, it is possible to show that each asset dynamics is equivalent to an individual instance of the hybrid Heston model. This allows, for example, to rely on existing techniques specifically devised for the Heston case, such as the numerical schemes for the solution of the Fokker-Planck PDE in the calibration of σ i (t, f ). More importantly, in the light of the above, we can interpret the WSLC model as an intuitive multi-asset extension of the Heston model able to introduce a high degree of flexibility in designing the dependence structure thanks to the underlying Wishart process.

Numerical results
This section is devoted to illustrate a realistic implementation of the WSLC model. We consider here a two-assets specification of the model intended to price derivatives on the pair of indices EuroStoxx50-DAX. We refer to Sect. 2.1 for the calibration of the two nonparametric local volatility surfaces (one for each asset, d = 2). Additionally, we consider the Wishart-based parameters in the third column of Table 2 in La Bua and Marazzina (2021), reported in the rightmost column of Table 4, i.e. we assume to deal with the parameters set obtained by requiring the Wishart process to be defined in the interior of the cone S + 2 (R), i.e. Condition (18). This setting is motivated by the fact that our goal is to construct a model able to produce realistic cross-asset correlations.
Having specified the LV and SV components, we can now address the task of calibrating the bivariate leverage function. To this extent, we estimate Eq. (35) by means of simulation-based approach (11) coupled with the GVA scheme described in La Bua and Marazzina (2021). This scheme can be extended to deal with the discretization of WSLC model trajectories by simply modifying Eq. (37) in La Bua and Marazzina (2021) intô for i = 1, . . . , d, where all other terms remain unchanged.
Despite the evidences found in La Bua and Marazzina (2021), indeed, it is our opinion that the GVA scheme should be preferred over the simpler TE scheme (also proposed in La Bua and Marazzina 2021) for calibration purposes. This is due to the fact that the GVA scheme allows for an accurate discretization of Wishart process. On the contrary, within the TE scheme we cannot bound the process (t) to remain in the cone of positive semidefinite matrices and diagonal elements can also become negative with nonzero probability. Even worse, the truncation mechanism could severely affect the estimate of (conditional) expected values in Eq. (35) and, ultimately, lead to wrong calibrated leverage functions. Having this in mind, we perform the calibration of leverage function by sampling 5 × 10 5 trajectories on a daily time grid. As in the single-asset case, we consider 50 equally weighted bins at each time step. We use the model so defined to price European call options written on each of the two indices, with maturity of 1 month and 3 years. As comparison, we price the basket of options in a pure LV setting and in the WASC model for which we consider the unconstrained parameters and those obtained when the condition β ≥ 3 is satisfied (see Table 4). In Table 5, we report the mean absolute error in volatility terms with respect to the market, while in Section B we exhibit the full set of numerical results.
The most evident result is that the WSLC model is able to overcome the limitations due to the request of a positive definite Wishart process. Not only the improvement with  Fig. 6 where the models induced implied volatilities are shown together with market ones. In particular, we want to stress that the only difference between the WSLC and the WASC model with the condition β ≥ 3 is the presence of the additional local adjustment factor in the former dynamics. As a consequence of Eq. (35), the hybrid model is in line with the LV one and reproduces satisfactorily the market implied volatility smiles, even when short times to maturity are considered. Overall, the new model succeeds in fitting the market evidences on single-asset markets and allows a sound modelling of the dependence structure.

Concluding remarks
In this article, we present a new class of hybrid models with the remarkable feature of embedding a matrix-defined dynamics for the stochastic evolution of variance factors. Numerical tests highlight that the new models effectively incorporate the advantages of LV approach in reproducing plain vanilla market evidences, still preserving the flexibility of Wishart-based pure SV models. In our framework, the contribution of the additional LV component is even more relevant if we take into account the fact that market calibrated parameters usually violate the condition for existence and uniqueness of (weak) solution to the Wishart SDE. We showed that the proposed SLV mixing presents an effective, probably unique, tool to deal with this problem within the class of Wishart-based pricing models. Additionally, WSLV model is found to be adequate even when we study the resulting volatility dynamics: it preserves the shape of forward implied volatility originated by the pure multidimensional SV model. This feature, in conjunction with the ability of properly manage the time structure of the skew, makes the WSLV a valid alternative to price forward implied volatility-dependent payoffs. Further, in the light of the evidences shown in La Marazzina (2019, 2021), it seems appropriate to consider WSLC as the most genuine and comprehensive framework proposed so far to extend the famous Heston model to the multi-asset case.
The resulting model, indeed, is able to generate a sophisticated dependence structure among assets and variance factors. More importantly, the new models succeed in eliminating any kind of mispricing of plain vanilla claims due to the stronger condition that we need to impose on model parameters in order to get reasonable cross-asset correlations dynamics. Calibration and pricing in the new framework are easily dealt thanks to the numerical techniques proposed in La Marazzina (2019, 2021). In particular, we extensively rely on the simulation schemes devised for WSV and WASC that allow to exploit the exact sampling of Wishart process in Ahdida and Alfonsi (2013).
In our opinion, the new models could turn out to be a comprehensive modelling framework to price a heterogeneous range of exotic derivatives consistently with European claims when a flexible multidimensional dynamics of variance factors is required.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Market indicates the implied volatility observed on 3 February 2016 for a given call option on DAX index. The error columns contain the difference between market implied volatility and the one obtained with the corresponding model Market indicates the implied volatility observed on 3 February 2016 for a given call option on DAX index. The error columns contain the difference between market implied volatility and the one obtained with the corresponding model