On the application of Wishart process to the pricing of equity derivatives: the multi-asset case

Given the inherent complexity of financial markets, a wide area of research in the field of mathematical finance is devoted to develop accurate models for the pricing of contingent claims. Focusing on the stochastic volatility approach (i.e. we assume to describe asset volatility as an additional stochastic process), it appears desirable to introduce reliable dynamics in order to take into account the presence of several assets involved in the definition of multi-asset payoffs. In this article we deal with the multi asset Wishart Affine Stochastic Correlation model, that makes use of Wishart process to describe the stochastic variance covariance matrix of assets return. The resulting parametrization turns out to be a genuine multi-asset extension of the Heston model: each asset is exactly described by a single instance of the Heston dynamics while the joint behaviour is enriched by cross-assets and cross-variances stochastic correlation, all wrapped in an affine modeling. In this framework, we propose a fast and accurate calibration procedure, and two Monte Carlo simulation schemes.

finance is fully motivated by the need to describe the multidimensional structure of asset variances (see La Bua and Marazzina (2019) and references therein), in this article our goal is to exploit the Wishart process in a multi assets framework. More precisely, we deal with the Wishart Affine Stochastic Correlation model (WASC), introduced in Da Fonseca et al. (2007) with the purpose of reproducing well-known multi-asset stylized facts in a tractable way. WASC model makes use of Wishart process to describe the stochastic variance covariance matrix of asset returns.
In our analysis we focus on the so-called Wishart processes introduced in Bru (1991) as a matrix generalization of square-root processes. A remarkable feature is that the analytical tractability is fully preserved since these processes belong to the class of affine processes. Given the strict connection with the well-known CIR processes, Wishart processes have been used to define multi-factor (Da Fonseca et al. 2008; La Bua and Marazzina 2019) and multi-asset (Da Fonseca et al. 2007) extensions of the classic Heston model, which is one of the most known and used models in finance, see for example (Goudenege et al. 2019;Yolcu-Okur et al. 2018).
Despite the analytical tractability (La Bua and Marazzina 2019), the implementation of Wishart-based models poses non-trivial challenges from a numerical point of view. In this article we extend the results presented in La Bua and Marazzina (2019) for the single asset Wishart Multidimensional Stochastic Volatility model (WMSV) to the WASC one. More precisely, we deal with the model calibration problem, presenting an innovative and efficient methodology to calibrate WASC parameters to market data. The algorithm exploits the close link existing between the Heston model (Heston 1993) and marginal WASC dynamics. Considering the single-asset case, in La Bua and Marazzina (2019) we show that, for an appropriate choice of parameters, both Heston (Heston 1993) and the Bi-Heston (Christoffersen et al. 2009) models may provide a reliable approximation of WMSV. In this article, for the multi-asset case, we extend this approximating technique making use of the distributional law of diagonal elements of Wishart process to connect the WASC calibration problem to the Heston one. We provide the analytical form of the gradient of calibration problem objective function with respect to Wishart-based parameters allowing for a further reduction in the computational burden.
Additionally, we propose model approximations that permit us to introduce two numerical schemes for Monte Carlo simulations. It is well known that a standard discretization (e.g. Euler scheme) is unfeasible in the Wishart-based framework, since we also need to take into account the evolution of non diagonal elements of the Wishart matrix (t) to determine the dependence structure and satisfy the positive semidefiniteness constraint for the Wishart process. Therefore, as a first algorithm, we propose an adapted version of the scalar full truncated Euler. Secondly, we extend the Gaussian variable approximation scheme presented in La Bua and Marazzina (2019) to the WASC model. Extensive numerical results to compare the two simulation schemes are provided.
The article is organized as follows. In Sect. 2 we present the Wishart process and its basic properties, while in Sect. 3 we deal with the WASC model. Finally the calibration procedure is described in Sect. 4, while the Monte Carlo simulation schemes are presented in Sect. 5. In both cases, numerical results are presented.

Definition of Wishart process and basic properties
In this section, we introduce the Wishart process. We refer to La Bua and Marazzina (2019) for details.
Definition 1 (Wishart process) Let W (t) be a d × d Brownian motion (i.e. a matrix of d × d independent scalar 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 stochastic differential equation (SDE): with , Q, M ∈ M d (R) (the set of real d × d square matrices).
As in La Bua and Marazzina (2019), in order to embed mean-reversion and stationarity, we consider matrix M to have only eigenvalues with negative real part. Moreover, we relate the deterministic part of the drift in (1), , to the expected long-term value of the process, ∞ , by the equation As shown in La Bua and Marazzina (2019), if we set d = 1 in (1), 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. Therefore Wishart processes can be considered as multidimensional extension of a scalar CIR process. Moreover, a Wishart process entails a non-trivial dependence structure among its elements, since we have where the notation [·, ·] refers as usual to the quadratic covariation of two stochastic processes, i j (t) indicates the element in the i-th row and j−th column of (t), and Q * = Q Q. Existence and uniqueness conditions of the solution of the Wishart process SDE are given by the following results.
Proposition 1 (Proposition 2.1 in La Bua and Marazzina (2019)) Let X (t) be a generic affine process with continuous trajectories defined in S + d (R) by the following SDE

) is a linear transformation. Such process admits a unique weak solution in
is the trace of a square matrix (i.e., the sum of the elements on the main diagonal).
If X (0) is in the set of real positive definite matrices S ++ d (R) and condition a) is replaced by the stronger requirement As stated in La Bua and Marazzina (2019), we can obtain the Wishart SDE (1) from (5) by setting D X = , C X = Q, and L [P 0 ] = M P 0 + P 0 M . Moreover, if we assume a more restrictive parametrization for the deterministic part of the drift conditions (a) and (c) of Proposition 1 are satisfied as soon as respectively, where the real positive parameter β plays the role of Feller's condition in the univariate case. Additionally if condition a) is not met the whole process is not well defined. For the rest of the paper we consider a Wishart process defined by (1) and (6) as usually done in financial literature. Notice that a significant constraint has thus to be imposed on parameter β. In the case d = 2, for example, we must require β ≥ 1. As shown in La Bua and Marazzina (2019), this condition is not usually met when we perform a straight calibration of Wishart-based pricing models to market prices of plain vanilla options, while it does not seem to be a real limitation, according to our knowledge, when a maximum likelihood or a moments estimation is considered (Alfonsi et al. 2016;Boloorforoosh et al. 2020;Da Fonseca et al. 2014;Gourieroux and Sufana 2010). Given our interest in exploiting the model for pricing purposes, we rely to the calibration on plain vanilla options, and therefore to the implied volatility surface, analyzing the impact of Conditions (7, 8) in our numerical results. In La Bua and Marazzina (2020) we show how this impact can be smoothed adding a local volatility component, i.e., exploiting a stochastic-local volatility hybrid model.

Distribution of Wishart process and related results
In this section we deal with the analogy between Wishart and CIR processes. In fact, as shown in La Bua and Marazzina (2019), exploiting the affine nature of the Wishart process, we have that its characteristic function is an exponential affine transformation of the initial state as shown in the following Proposition: Proposition 2 (Proposition 2.2 in La Bua and Marazzina (2019)) Let be a real symmetric d × d matrix, t ≥ 0 and T − t = τ > 0. The (conditional) characteristic function of the Wishart process defined by (1) and (6) is where ι is the imaginary unit (i.e. ι = √ −1) and matrix A ( , τ ) and scalar function b ( , τ ) are such that The additional matrix functions appearing in the above equations are given by As a consequence of the analytical tractability of Wishart process and of the knowledge of its characteristic function, we are able to present an additional result regarding the distribution of diagonal elements of the Wishart process (a similar result for the trace of the Wishart process has been obtained in La Bua and Marazzina 2019, Corollary 2.3).
Corollary 1 (Distribution of elements on the main diagonal of Wishart process) Let i (t) = ii (t) be the i-th element on the main diagonal of (t) and F χ 2 (x; ν, δ) the cumulative distribution function of a non-central chi-square random variable with ν degrees of freedom and non-centrality parameter δ. Then, for a fixed T > t, we have: where δ i = γ i /ϑ i with (τ ) = (γ i j ) 1≤i, j≤d and (τ ) = (ϑ i j ) 1≤i, j≤d . For ease of notation, we also set γ i = γ ii and ϑ i = ϑ ii .
Proof We use the fact that exp Tr log(G) = det [G] for any matrix invertible matrix 1 G to write Let λ be a real variable and e d i = (1 k= =i ) 1≤k, ≤d , then by setting i = λe d i , the characteristic function of i (T ) is from which (10) follows by the definition of non-central chi-square distribution.
An alternative way to obtain (10) is to recognize that (11) is the characteristic function associated to the non-central Wishart distribution with scale matrix (τ ) and non-centrality matrix (τ ) and apply results in Kourouklis and Moschopoulos (1985). For the sake of completeness, we point out that an analogous claim is shown in Da Fonseca and Grasselli (2011). Our formulation, however, gives a direct interpretation of parameters involved in the distribution of i (T ) in terms of matrices describing the Wishart process. As we will see, this turns out to be particularly useful for computational purposes.
An important consequence of (10) is that we can define an exact mapping between i (T ) and a CIR process: Proposition 3 (CIR process mapping i (T )) Let v(t) be a CIR process defined by (3). For a fixed T > t, it holds that (conditionally on v(t) and (t) respectively) v(T ) and where γ i and ϑ i have been introduced in Corollary 1.
Proof The correspondence of the distributions relies on the properties of the CIR process. We refer to Cox et al. (1985) for details.

The wishart affine stochastic correlation model
With the purpose of reproducing well-known multi-asset stylized facts in a tractable way, in Da Fonseca et al. (2007) the authors introduce the Wishart Affine Stochastic Correlation model (WASC) that 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: where diag [·] is the operator that transforms a d-dimensional column vector into a with z(t) another d-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 (18) 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 (17) allows for stochastic correlation among asset returns in a tractable framework where each asset is enriched with a stochastic volatility behavior consistent with the effects observed on plain vanilla markets. Let the i-th forward asset price at time t be denoted by f i (t) = f i (0)e y i (t) , t ≥ 0. The peculiarities of the model can be fully appreciated by referring to the individual, or scalar, dynamics of asset returns y i : where S(t) = S(t) = s i j 1≤i, j≤d is the unique positive semi-definite square root of (t). We also use the notation i (t) = ii (t) to denote the i-th diagonal element of Wishart process. By straightforward computations, from (19) we can compute the quadratic covariation of two given assets: that highlights the role of Wishart process, used to describe the stochastic evolution of the asset returns variance covariance matrix. Furthermore, we can explicitly define the cross-asset correlation matrix C y (t) as By exploiting the properties of Wishart process, it can be shown that C y (t) is a welldefined correlation matrix (i.e. C y (t) is positive semi-definite and each ρ i j (t) ∈ [−1, 1]) as soon as condition (7) is satisfied and provided that i j (t) = 0 for we use the following theorem that applies for Wishart distributed matrices: Theorem 1 (Theorem 2.4.2 in Kollo and von Rosen (2006)

symmetric matrix that follows a non-central Wishart distribution with degrees of freedom β, scale and non-centrality matrix . Then for
If we set X W = (t) and B = e d i , e d j for some admissible i and j (with e d i the i-th element of the standard basis of R d ), the previous theorem shows that the resulting is a well-defined 2 × 2 Wishart process and By combining (19) and (20) with Proposition 3, and fixing a time horizon T , we can represent the (T -specific) WASC dynamics as the following 2d system of scalar where the parameters κ i , θ i and η i are given, respectively, in (14), (15) and (16). Here the correlation structure among Brownian motions w = w 1 , w 2 , ..., w d , w y 1 , w y 2 , ..., w y d is described by means of the stochastic block matrix where the submatrices, other than C y already introduced in (21), will be described in the following. Interestingly, from (22) and (23) we have that, for any T > 0, the scalar dynamics of each asset is consistent with a standard Heston model driven by the i-th diagonal element of (t). This assures that the behaviour induced by WASC in terms of reconstructed implied volatility surfaces is in line with the documented findings of traditional one factor stochastic volatility models. Consistently with Heston model, the asset specific returns-volatility correlation is constant: following Da Fonseca et al. (2007), we have where R i is the matrix with r on the i-th row and zero elsewhere. We can even generalize the previous result by explicitly deriving the correlation between the i-th log-asset and a generic diagonal element of (t): we, indeed, have that (as shown in Da Fonseca et al. (2007)) the covariation between asset returns and volatility terms is given by By combining (26) and (20) and exploiting (4), we get the generic element of matrix C y (t) as In the last equality of (27), we introduce a new representation of such correlations that can be seen, quite fascinatingly, as the product of the (constant) proper, or scalar, j-th asset-volatility correlation and the cross-asset correlation between i-th and j-th assets. This result highlights the peculiar dependence structure inherent in the WASC model and could help in gaining more insights on parameters impact on correlation surfaces in the spirit of the study carried out in Da Fonseca et al. (2007). Further, using (4), it follows that the elements of C (t) have the form From (21), (27) and (28), we have that the stochastic evolution of (24) is fully described by processes ρ i j (t). Let us consider, for example, the case d = 2: matrix C(t) then reads as . This representation highlights the peculiar dependence structure induced in the WASC model and could provide useful insights on the role of Q and r in determining the relation among state variables.

A restricted version of the model
In this section we consider a restricted, more intuitive, specification of WASC model: we assume matrix M to be diagonal and with negative entries. This setting leads to a very interesting dynamics for the diagonal elements of Wishart process. From Proposition 3, if M is diagonal, direct computation gives for all T > 0. An immediate consequence is that the asset instantaneous variances are now described by time-independent CIR processes, in the sense that the parameters involved are no longer function of the time horizon considered. In our opinion the resulting parametrization turns out to be the most genuine multiasset extension of the Heston model: each asset is exactly described by a single instance of the Heston dynamics while the joint behaviour is enriched by cross-assets and crossvariances stochastic correlation, all wrapped in an affine framework. As far as we know, there are no alternative settings that can reach a comparable degree of flexibility. The exact Heston representation of asset dynamics also helps in understanding the role and the impact of WASC parameters, that in the general formulation appear somehow unclear. In particular, it is worthwhile to point out that the pricing of single-asset European claims is only affected by the corresponding diagonal element of 0 . To see this, it suffices to notice that for λ = λ i e d i , the matrix A y (τ ) in (31) has a non-null i-th diagonal element and zeros elsewhere. This peculiarity is in line with the asymptotic analysis provided in Da Fonseca and Grasselli (2011) where the i-th implied volatility approximation for short time to maturity is found to be 2 denoting the log-forward moneyness. Consequently, the off-diagonal entries of 0 can be used to match multi-asset stylized facts without compromising the shape of individual volatility surfaces. This represents an additional degree of freedom that in the general WASC model we would not have. A possible calibration strategy could be to set the off-diagonal entries of 0 in order to match a predefined initial cross-asset correlation matrix. Alternatively, provided that liquid multi-asset derivatives are traded, we could try to fit the implied correlation market evidences. In order to further develop this point, we now study the impact of 12 on the price of Best-Of put options which payoff is S i (t) being the asset price at time t, t ≥ 0. Let π W A be the following set of WASC parameters: that are meant to describe realistic market scenarios and allowing for a well defined Wishart process. Figure 1 shows the implied correlation profiles corresponding to different values of 12 , where the implied correlation is defined to be the value of parameter ρ such that the WASC price equals the one obtained in a two-assets Black-Scholes setting, i.e.
Here σ W A imp,i (K , T ) is the Black-Scholes implied volatility corresponding to the option written on the i-th asset with strike K and maturity T whose price is computed with WASC model. By exploiting the affinity of the model, Best-Of put options are priced by numerically computing a bi-dimensional inverse Fourier transform. We refer the reader to Da Fonseca et al. (2007) where this pricing methodology is developed for Best-Of contracts. From the numerical results, it is evident that the off-diagonal Wishart element plays a significant role in modelling the implied correlation skew: we observe an increase in implied correlation levels for higher values of 1,2 . This, in turn, induces an increase in option prices consistently with the fact that Best-Of put options are long correlation products that benefit from lower assets returns dispersion.

WASC Characteristic function
The chosen correlation structure (18) assures the affinity of WASC model. This means, once more, that we can express the (joint) characteristic function of the asset returns vector y(T ) as an exponential affine transformation of state variables y(t) and (t) as recalled in the following Proposition: Proposition 4 (Joint characteristic function of log-prices in WASC) Let the logforward prices vector y(t) be described by (19) with the deterministic matrix A y (τ ) and the scalar function b y (τ ) given by .
Thanks to Proposition 4, we are able to price both plain vanilla and multi-asset options (if transform-based techniques are applicable 3 ) in a comprehensive framework. In particular, we price options on the i-th asset as a basket option with degenerate weights vector e d i , such that λ = λ i e d i . Furthermore, the knowledge of the joint characteristic function of asset returns vector allows to make use of bounds techniques as those developed in Caldana et al. (2016) for basket options. Despite the analytical tractability, several numerical issues arise when we try to calibrate WASC model to market data by exploiting (31). Not only, indeed, we have to evaluate functions of matrix argument for each computation of the characteristic function (as in the WMSV case), but, even worse, we are required to perform d different plain vanilla pricing (one for each asset) for a single parameters set. This is due to the lack of liquid multi-asset derivatives that force us to calibrate model parameters to the individual market implied volatility surfaces. As reported in Da Fonseca and Grasselli (2011), such a naive algorithm can take up to 15 minutes in the simplest case d = 2. This is, clearly, not feasible for real market applications. Therefore in the next section we propose an accurate and fast calibration procedure.

A new calibration procedure
In this section we present an innovative and efficient methodology to calibrate WASC parameters that exploits the close link existing between Heston model and marginal WASC dynamics. The proposed algorithm is firstly tested in a simplified framework and then applied to market data. The results obtained also highlight the impact of parameter β on model accuracy in reproducing market volatility smiles.
Let us consider a WASC parameters set π W A (with cardinality N W A ) and fix a maturity T . For the generic i-th asset described by (19) (13)-(16) Then, the Jacobian matrix ofr i,T (π W A ) (the residuals vector composed of options with maturity T written on the i-th asset) can be written as where the second matrix in the right-hand side of (33) is known thanks to Cui et al.  (2019), the gradient of the calibration problem objective function is given by ∇ f obj = J W Ar (π W A ).
The calibration algorithm so defined avoids the computation of WASC characteristic function and significantly reduces the issues due to the possible presence of multiple minima. Given that we rely on the law identity in Proposition 3 rather than on some approximation, the routine does not require any further step.

A simplified calibration exercise
The accuracy of the proposed algorithm is illustrated by considering the following numerical experiment: let us suppose that a fictitious two-assets market is perfectly described by the WASC parameters reported in the first column of Table 1. Even if simplified, the data outline realistic market environments: they represent the calibrated parameters set (truncated at the first significant decimal digit) found in Da Fonseca and Grasselli (2011) for the couple of indices EuroStoxx50-DAX. We construct a full implied volatility surface for each asset assuming to have options with maturities T = [0.25, 0.5, 1, 3] and 41 equally spaced strikes ranging from 0.5 to 1.5 (initial asset values are set for simplicity equal to 1). Each of the resulting surfaces consists of 164 options. The goal is to implement the proposed algorithm in order to find a suitable parameters set that reproduces the supposed market data. Hopefully, we expect the calibrated parameters to be reasonably close to the original ones (accuracy) and to experience a limited dependency on the initial guess (robustness). For this test we set the starting values of the optimization routine as shown in the second column of Table 1. The choice is meant to assess the robustness of the algorithm in the case in which the initial guess is very far from the optimal set. Indeed, not only the discrepancy is mixed -some values are overestimated, others underestimated -but the distance between initial guess and optimal values is substantial: the smallest gap, defined as percentage difference, is equal to 35.29%. The mistaken initialization of the problem and the high dimensionality of the parameters space make the calibration task more challenging and could potentially lead to suboptimal outcomes. Notwithstanding, the proposed algorithm is able to produce results very close to the original values: the norm of the errors between true prices and calibrated ones is 2.2069 × 10 −7 . Most remarkably, the procedure takes only 3.56 seconds using Matlab Mex files, on a laptop PC with an Intel Core i7 CPU and 8 GB RAM. By considering parallelization and porting to more efficient languages we can obtain a further speedup. It is worthwhile also noting that in realistic applications, the calibration problem is somehow facilitated thanks to the availability of previous optimal sets that act as efficient guesses. In the lights of all these evidences, we believe that the proposed methodology represents a highly efficient tool for the calibration of WASC model. This is particularly true if we intend to increase the number of assets involved with the subsequent growth of dimensionality.

Calibration to market data
We now want to validate the procedure with realistic market data. Despite the general applicability of the algorithm, we focus our attention on the restricted specification of the model introduced in Sect. 3.1. With this in mind, we select a basket of market quoted instruments composed of 201 European call options written on EuroStoxx50 index and 182 on DAX index. The set of derivatives on the DAX is the same set  (2019). We further set, for simplicity, interest rates and dividends to zero. Thanks to the efficiency of the new calibration algorithm, we are able to calibrate model parameters in less than 3 seconds. The outputs of the optimization routine are shown in the leftmost column of Table 2. Given that, as illustrated above, the off-diagonal element of 0 does not impact the pricing of univariate call options, we set its value such that the initial correlation among the two indices equal the one-year historical one (that is found to be 0.9715). 4 The most interesting result is that β is lower than 1. This is coherent with the evidences in Da Fonseca and Grasselli (2011) where similar results are found. Figure 2 shows the calibrated model implied volatility skews for the two indices with respect to maturities of one month, one year and three years. The model succeeds in reproducing the shape of market volatility surfaces but the mispricing is not negligible for short term far-from-the-money options. This is particularly true for the EuroStoxx50 index as highlighted from the fact that the error in volatility terms is roughly 3 times higher than the error made for the DAX.
Additionally, we can compare the evidences from the calibration of the WASC model against the Wishart Multidimensional Stochastic Volatility model (WMSV, La Bua and Marazzina 2019), as well as the Heston (1993) and Christoffersen et al. (2009) models, calibrated to the same basket of DAX options. Table 3 shows the calibrated initial variance of asset returns along with the Mean Squared Error with respect to both price and implied volatility for the four models. Consistently, the estimates of initial variance are in strict agreement: all the models agree on the initial volatility. Moreover, w.r.t. the accuracy of the calibration, the two multi-factor models, i.e., the WMSV and the Bi-Heston, tend to perform quite similarly (although errors for WMSV are slightly smaller) and substantially outperform the simpler Heston and  The Initial Variance row specifies the initial variance of asset returns in the corresponding model WASC dynamics. In particular, by comparing the error of WMSV and WASC models, the outperformance of the former is clearly evident. This is not surprisingly since we contrast a multi-factor volatility setting (WMSV) with the WASC single-asset dynamics that, as developed in Sect. 3, is equivalent to 1-factor parametrization. It is important to remark, however, that the models are meant to address rather different tasks (i.e. single-asset and multi-asset modelling). Moving back to the multi-asset calibration in Table 2, some fix is required in order to enforce the existence and uniqueness condition for matrix process (t). First of -   Table 2. Even if the loss in accuracy seems to be somehow limited (the error measure are just slightly higher than in the unconstrained setting), a relevant issue arises: in Fig. 3 we report simulated trajectories of WASC cross-asset correlation obtained with the resulting parameters set. The fact that β satisfies condition (7) effectively ensures ρ 12 (t) to lie in the range [−1, 1]. However, the fact that the parameter is just slightly above the threshold (β = 1) makes the boundary of S + 2 (R) very likely to be attained. Very often, then the absolute value of correlation is stuck at 1. We can also experience sudden changes in correlation from +1 to −1 (or viceversa) in a very restricted time frame (even on a daily basis). In order to study the dependence of the observed phenomenon on the initial value of correlation, in the rightmost panel of Fig. 3 we also consider the case 12 = 0 that produces a similar erratic correlation dynamics.
Given the intent to apply WASC model to describe the joint behaviour of asset prices, this represents a major issue. To tackle the problem, we decide to enforce the positive definiteness condition for (t), given by (8), that in our setting equals to set β ≥ 3. Calibrated parameters are collected in the rightmost column of Table 2. The corresponding cross-asset correlation dynamics is depicted in Fig. 4: the trajectories are now much more meaningful. Further, as a consequence of the fact that (t) is defined on the interior of S + 2 (R), ρ 1,2 (t) is bounded in (−1, +1). Nonetheless, the stronger condition enforced has a severe impact on the ability of the model to reproduce singleasset market evidences. Reconstructed volatility skews are shown in Fig. 5. Significant discrepancies now emerge for far-from-the-money options. In particular, in the very short-end of the volatility term structure the error with respect to market volatilities can be as high as 11.87% (in-the-money options on EuroStoxx50) and 10.80% (out-ofthe-money options on DAX). Disappointingly, we face a non trivial trade-off between plain vanilla pricing accuracy and realistic modellization of cross-asset correlation. A possible solution to mitigate the problem could be to set β equal to some value in the range (1, 3]. This alternative, however, would require to couple the plain vanilla analysis with adequate market evidences on multi-asset derivatives.

Simulation schemes for the WASC
This section is devoted to present simulation algorithms specifically devised for WASC model. As far as we know, indeed, there are no previous attempts in literature to deal with the discretization of prices trajectories (17). In particular, our task is to develop an efficient, yet accurate, scheme to discretize the system of SDEs (22)-(23). It is evident that a standard discretization (e.g. via Euler scheme) is unfeasible, since we also need to take into account the evolution of non diagonal elements of (t) to determine the dependence structure and satisfy the positive semi-definiteness constraint for the Wishart process.  As a first algorithm, we implement an adapted version of the scalar full truncated Euler (TE) scheme that reads as with z and W , respectively, d-dimensional vector and square matrix of independent standard gaussian random variables. Here Vec [·] is the operator that extracts the elements on the main diagonal of a square matrix into a column vector, while + is the positive part of matrix . As a second algorithm, we deal with the Wishart process sampling scheme developed in Ahdida and Alfonsi (2013): that is, we consider as given an entire discretized path of (t) over the time grid 0 = t 0 < t 1 < ... < t M T = T with time step . In this way we are only left with the problem of sampling the log-prices trajectories. The most challenging task here is to embed the correlation structure (24) in the discretization of (22). In standard cases, we would compute the Cholesky decomposition of matrix C(t) such that where L C (t) = i, j (t) 1≤i≤d,1≤ j≤i is the lower triangular matrix that satisfies is a vector of independent Brownian motions. By exploiting (34), we can rewrite (22) as that can be discretized by generating 2d independent gaussian random variables. Unfortunately, in our setting this is not readily doable as a consequence of the mechanics of Wishart sampling algorithm. In other words, we do not have a direct "access" to the discretized paths of w 1 , w 2 ,..., w d . The simple idea underlying the new simulation scheme is to exploit the auxiliary scalar dynamics of i (t) to get an approximation of w i (t + ) − w i (t). Letˆ (t) andˆ (t + ) be the realizations of the trajectory of Wishart process for two adjacent points on the time grid computed by means of the exact scheme in Ahdida and Alfonsi (2013). The discretized version of (1), that readŝ can now be used to approximate the gaussian variable w i . Let w i be the result of (36), for a sufficiently small time interval, w = w 1 , w 2 , ..., w d represents an approximation of a vector of gaussian variables with correlation matrixĈ (t) (i.e. the realization at time t of matrix C ). Further, letL (t) be the lower triangular matrix obtained from the Cholesky decomposition ofĈ (t), then is composed of d approximated independent gaussian random variables. By sampling an additional random vectorŵ * y from N (0 d , I d ) (the d-variate gaussian distribution) and settingŵ * = w * ŵ * y we can finally approximate (35) aŝ whereL C (t) results from the factorization ofĈ(t). IfĈ(t) turns out not to be positive definite, we take its positive part,Ĉ + (t), (defined as the matrix obtained from the spectral decomposition ofĈ(t) with negative eigenvalues replaced by zeros) and apply the extended Cholesky decomposition described in Golub and Van Loan (2012). The complete algorithm is exhibited in Algorithm 1, and we refer to it as Gaussian Variables Approximation (GVA) scheme for WASC model. Notice that this algorithm can be considered as an extension of the GVA algorithm for the WMSV case developed in La Bua and Marazzina (2019).

Numerical results
Even if these new schemes would apply to the general specification of the model, here we suppose to deal with the reduced model presented in Sect. 3.1 (i.e. we consider matrix M to be diagonal). In particular, we develop an extensive numerical investigation based on the parameters set calibrated to market data enforcing the condition β ≥ 3 and shown in the rightmost column of Table 2. Considering 5 × 10 5 simulation paths, we price European call options written on any of the 2 assets with maturity T = 1 and moneyness in the range {70%, 100%, 130%}. For the sake of simplicity, we assume interest rates and dividends equal to zero and f 1 (0) = f 2 (0) = 100. The asterisk in the following tables means that the corresponding reference value lies outside of the 95% confidence interval.
In the context of plain vanilla options, results in Tables 4, 5, 6, 7 show that both schemes allow for accurate price estimates as the size of the time step is sufficiently small. More in detail, the TE scheme is found to outperform the GVA scheme when the time grid is coarse (10, 20 and 50 steps per year), while for smaller mesh widths the two approaches tend to perform similarly. With this setting and taking into account option prices for both assets, the absolute mean percentage error is respectively equal to 0.232% for the GVA scheme and to 0.275% for the TE scheme. It is worthwhile to remark, though, that the GVA scheme systematically requires a finer time discretization to produce reliable estimates (true prices lying in the 95% confidence interval) compared to the simpler TE scheme. This is due to the fact that the approximation exploited in (36) seems to be adequate only for very small time intervals. From a computational point of view, the TE scheme greatly outperforms the GVA scheme with the latter that results 110% − 130% slower than the former as exhibited in Table 8. When implementing the GVA scheme, indeed, at each time step we are asked to perform the Cholesky factorization of the 2d-dimensional matrix C(t) with a considerable increase in the computational burden.    We now compare the two simulation schemes pricing an Asian option with payoff in T = 1 equal to with M = 1 12 , i.e., we are considering an Asian option with monthly monitoring, and K = 200 (all the other parameters as above). In Table 9 we compare the two simulation schemes: results confirm that the TE scheme outperforms the GVA when the time grid is coarse.
Nonetheless, the GVA scheme proposed embeds the inherent advantage to implement the exact sampling of Wishart process thanks to the algorithm in Ahdida and Alfonsi (2013). This feature is of great importance in all the cases in which we need to estimate the (conditional) moments of the distribution of the elements of (t). The ability to consistently deal with the discretization of Wishart process is, indeed, the main reason that led us to develop the new scheme. To better highlight this advantage, in Table 10 we exploit the two simulation algorithms to compute the absolute value of the real part of the joint characteristic function setting X(T ) = log( f 1 (T )) log( f 2 (T )) , X = 0 a , V = 1 0.5 0.5 1 , varying the parameter a. Notice that a closed form solution for this characteristic function is provided in Da Fonseca and Grasselli (2011). For large values of a, e.g., a = 1, the TE scheme outperforms the GVA one, exactly as in the option pricing cases above described. However, the opposite happens for small values of a, i.e., when the contribution of the simulation of the Wishart process is more important w.r.t. the logasset value. Moreover, when a = 0, i.e., we only need to simulate the Wishart process (T ), the GVA scheme performance is independent on the number of the time steps, since it is an exact simulation scheme, as shown in Ahdida and Alfonsi (2013).

Concluding remarks
The matrix structure of Wishart-based stochastic volatility models provides a remarkable degree of flexibility in describing the evolution of asset(s) volatility. Realistic implementations, though, require the development of specific numerical techniques in order to deal with the inherent level of complexity. In this article we have shown, leveraging on a thorough analysis of distributional properties of Wishart process, some possible solutions intended to make this class of model more suitable for real market applications. Accordingly, we hope that our contribution will increase the interest of researchers and practitioners towards matrix-variate stochastic volatility dynamics.
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/.
The unspecified elements of J H −W A are equal to zero.