Variance and Dimension Reduction Monte Carlo Method for Pricing European Multi-Asset Options with Stochastic Volatilities

Pricing multi-asset options has always been one of the key problems in financial engineering because of their high dimensionality and the low convergence rates of pricing algorithms. This paper studies a method to accelerate Monte Carlo (MC) simulations for pricing multi-asset options with stochastic volatilities. First, a conditional Monte Carlo (CMC) pricing formula is constructed to reduce the dimension and variance of the MC simulation. Then, an efficient martingale control variate (CV), based on the martingale representation theorem, is designed by selecting volatility parameters in the approximated option price for further variance reduction. Numerical tests illustrated the sensitivity of the CMC method to correlation coefficients and the effectiveness and robustness of our martingale CV method. The idea in this paper is also applicable for the valuation of other derivatives with stochastic volatility.


Introduction
In the last 40 years, financial derivatives have become increasingly important in finance.They are actively traded on many exchanges throughout the world, and are entered into by financial institutions, fund managers, and corporate treasurers in the over-the-counter market.They are especially important for market anticipants because they can be used to transfer a wide range of risks in the economy from one entity to another.Efficient use of financial derivatives can certainly promote financial and social sustainability.For instance, there are many different types of energy and agricultural commodity derivatives that are designed and used to contest weather and market risks, and to protect the benefit and reduce the potential loss of anticipants.Another example is that the real options approach is very popular in valuing the real estate sustainable investment.Conversely, inappropriate use of derivatives may cause great, even global, disasters, for example, the credit crisis that started in 2007.As pointed out by Hull [1], "we have now reached the stage where those who work in finance, and many who work outside finance, need to understand how derivatives work, how they are used, and how they are priced." The accurate and fast pricing of financial derivatives is one of the most important things in financial engineering since many of the problems in economics and finance eventually turn into the pricing of financial derivatives.For example, Kim et al. [2] decided the optimal investment timing using rainbow options valuation for economic sustainability appraisement.Yoo et al. [3] determined an optimum combination of financial models including options to achieve a sustainable profit for overseas investment projects.The pioneering work of Black and Scholes [4] and Merton [5] lay the foundations for option pricing models.It is well known that the stochastic volatility model can be used to generalize the constant volatility assumption in the Black-Scholes model to capture the character of empirical observations from financial markets, such as the observed volatility smile and the leptokurtic features of the asset return distribution [6,7].Stochastic volatility models describe volatility behavior with another stochastic differential equation.There are many studies on stochastic volatility models, such as those of Hull and White [8], Scott [9], Stein and Stein [10], Ball and Roma [11], Heston [12], Schöbel and Zhu [13], and Hagan et al. [14].In addition to these one-factor stochastic volatility models, Fouque et al. [15][16][17][18] proposed a multi-factor mean-reverting stochastic volatility model.A comprehensive treatment of stochastic volatility models can be found in Reference [19].
Multi-asset options refer to a wide variety of contingent claims whose payoff depends on the overall performance of more than one underlying asset.Usually, they can be grouped into three categories: rainbow options, basket options, and quanto options.The prices of rainbow options rely on price changes of underlying assets, such as exchange options, outperformance options, spread options, chooser options, max-call options, and their variations.Basket options prices are always determined by the average price of underlying assets, while the value of a quanto option depends on the performance of domestic and foreign underlying assets.Jiang [20] introduced the concepts and constructed the pricing models of multi-asset options in detail, where the volatilities were constant.The pricing problem of multi-asset options pricing is essentially equivalent to a high dimensional integral.It is challenging to compute such a high dimensional integral, especially when the number of underlying assets is large, or stochastic volatilities are considered in the model.
There are mainly three pricing methods for multi-asset options: the analytic approximation method, the fast Fourier transformation (FFT) method, and the MC simulation method.The analytic approximation method typically constructs an approximate pricing model for the original problem that results in a closed form solution, and this method is always elegantly designed to the original problem.Several studies focus on this approach, for instance, those by Turnbull and Wakeman [21], Curran [22], Milevsky and Posner [23], Ju [24], Zhou and Wang [25], Alexander and Venkatramanan [26], Datey et al. [27], Brigo et al. [28], Borovkova et al. [29], Deelstra et al. [30,31], and Li et al. [32].The main disadvantage of the analytic approximation method is that the size of the error is unknown and there is no way to systematically reduce it.The FFT method, which was prosed by Carr and Madan [33], has successfully been used in option pricing problems with a low dimension because of its high efficiency (see Carr and Wu [34], Heston [12], Grzelak et al. [35][36][37], and He and Zhu [38]).However, the FFT method depends on the availability of a characteristic function (usually in an affine framework), which is not always promised in a general stochastic volatility model.It is also difficult to apply the FFT method to high dimensional problems due to their dimensionality.Thus, for higher dimensional options, the most practical method seems to be MC simulations.Kim et al. [2] and Yoo et al. [3] also used MC simulations to price the embedded option prices in valuation real investment projects since the high dimension of problem.MC uses the sample mean as an estimator for the expectation of a random variable.Its speed of convergence is not influenced by the dimension of the problem.In addition, it allows for a simple error bound, given by the central limit theorem.
The major drawback of an MC simulation is that its convergence rate is quite slow, that is, O(m −1/2 ), where m is the number of samples in MC simulation.As a result, often the main challenge in developing an efficient MC method is to find an effective variance reduction technique.There are a lot of studies about how to improve the efficiency of an MC simulation, and we refer the reader to Glasserman [39] and relevant references therein for a detailed discussion on various variance reduction techniques.Kemna and Vorst [40] presented one of the classical works on accelerating MC simulations.They used the geometric average option as a CV to price the arithmetic average option, which proved to be very successful.For a European multi-asset option pricing problem, Barraquand [41] proposed a "quadratic resampling" method by matching the moments of the underlying assets to reduce the variance of the MC simulation.Pellizari [42] designed a CV method called mean Monte Carlo to gain variance reduction of an MC simulation.The key of their success was that a Black-Scholes formula could be obtained when all underlying assets except for one were replaced by their mean.Borogovac and Vakili [43] proposed a "database Monte Carlo" CV method that avoids computing the expectation of CV, but the database, constructed in advance, requires huge calculations.Dingeç and Hörmann [44] exploited the property that the geometric average price was larger than the arithmetic average price to construct a CV by conditioning the payoff on the assumption that the geometric average of all prices was larger than the strike price.The expectation of their CV was computed by numerical methods, and their numerical tests for Asian options and basket options showed a great accelerating effect on the MC simulations.Liang et al. [45] designed a CV for European multi-asset options based on principal component analysis.Sun and Xu [46] used the CMC method with the importance sampling technique to accelerate MC simulations for basket options.There are some other approaches to speed up MC simulations, such as the quasi-Monte Carlo method [39,[47][48][49][50][51][52], and parallelized implementations of MC simulations on CPUs/GPUs [53][54][55][56][57][58].
However, there is little research on variance reduction of MC simulations in pricing multi-asset options with stochastic volatilities.Du et al. [59] proposed a variance reduction method in multi-asset options under stochastic volatility models by matching the moments of the volatilities.Although their method shows great variance reduction of MC simulations, there are some restrictions to it.(1) All underlying assets are assumed to be driven by one stochastic volatility factor, which is not reasonable in practice.A more reasonable model is to assume that each underlying asset is driven by its own stochastic volatility factor (see Antonelli et al. [60], Shiraya and Takahashi [61], and Park et al. [62]).
(2) Their moment matching technique greatly depends on the Hull-White stochastic volatility model, and is not applicable to general stochastic volatility models.(3) They only conducted numerical tests for options with two assets, which is not general enough for most multi-asset options.
In this paper, we aim to develop an efficient dimension and variance reduction method for MC simulations in pricing European multi-asset options with general stochastic volatilities.In our pricing framework, the underlying asset is assumed to be driven by its own stochastic volatility process, and full correlations between factors are allowed.The stochastic volatility model, which could be the Hull-White [8] or Heston [12] models, is quite general, such that our pricing model has a wide range of applicability.Our dimension and variance reduction method is built on the idea developed by Liang and Xu [63], who designed a CMC simulation with martingale CV to price single-asset European options with stochastic volatility.Our main contributions are: (1) A CMC pricing framework is deduced for European multi-asset options with general stochastic volatility models, which results in dimension and variance reduction.(2) A martingale CV based on a martingale representation theorem is combined with the CMC to obtain further variance reduction of the MC simulations.(3) The algorithm was tested on typical multi-asset options, such as exchange options, basket options (which can be more than two assets), and quanto options, showing the broad applicability and high efficiency of our method.
The paper is organized as follows.In Section 2, we introduce the pricing model of multi-asset options with stochastic volatilities.In Section 3, we deduce the CMC pricing framework, prove the martingale presentation theorem, and construct the martingale CV in detail.We present numerical tests and their results in Section 4, to evaluate the efficiency of our proposed method.Finally, we conclude the paper in Section 5.

Pricing Model
In this section, we give the pricing model of multi-asset options with stochastic volatilities.Specifically, in a risk-neutral world, let S i (t) be the price of the ith underlying asset (i = 1, 2, • • • , n) at time t, which we assume obeys the following stochastic differential equations: where r is the constant risk-free interest rate and δ i is the continuous dividend rate.Y i is the stochastic variance, and the functions f i (y), µ i (t, y), and g i (t, y) determine the specific volatility model, which can be quite general.dW i (t), and dZ i (t) are standard Brownian noise terms, and the covariance between them is captured as follows: where the correlation coefficients ρ ij , ρ i are constant.Equation (3) indicates that the underlying assets are correlated.Equations ( 4) and ( 5) indicate that any underlying asset is driven by only one stochastic variance factor and is not directly affected by the other stochastic variance factors.Equation (6) indicates that the random stochastic variance factors are mutually independent, but this assumption could be relaxed allowing for correlated random processes.Several popular stochastic volatility models are collected in Table 1.
Hagan et al. [14] = 0 y 0 σy Notes: µ is the drift of Hull-White stochastic volatility model.σ is the volatility of stochastic volatility.a is the rate of mean reversion and θ is the long-term mean of stochastic volatility.All parameters µ, a, θ, σ here are constants.The functions µ i (t, y), g i (t, y) here have no explicit dependence of time t.
In the following, we introduce our notations for convenience.The underlying asset vector is S(t) = (S 1 (t), • • • , S n (t)) , and the stochastic variance vector is , where represents the transpose of a vector or matrix.Additionally, the correlation matrix is given by Γ = (ρ ij ).
Suppose the payoff function of the European multi-asset option at maturity T is given by: Denote V(t, s, y) as the value of a European multi-asset option with stochastic volatilities at time t; then, by the no-arbitrage pricing principle we obtain: where E[•] is the expectation in a risk-neutral world.Given the initial asset price S(0) and initial variance Y(0), the European option price at the initial time is actually: MC simulation can be used to compute the option price based on Equation (9) (please see Glasserman [39]).Suppose the number of samples in MC simulation is m.Firstly, for the jth sample, we need to simulate the processes of the Brownian motions W Then, we average the samples of discounted payoff and use the sample mean V = 1 m ∑ m j=1 V (j) as an estimation of the option price.The law of large numbers guarantees the convergence of MC simulation.The central limit theorem guarantees that the standard error-the standard deviation of sample mean V-from MC simulation has a form of Std = var(V (j) ) √ m .The standard error can be used to measure how far the sample mean is likely to be from the option price or to make confidence intervals of the option price, for instance, a 95% confidence interval V ± 1.96Std.It also indicates that the MC simulation has a convergence rate as O(m −1/2 ), which is rather low.Thus, in the next section, using a similar idea as in Liang and Xu [63], we propose our efficient CMC simulation framework with martingale CV for this option pricing problem.
If the stochastic volatility f i (Y i (t)) in Equation ( 1) is replaced by a constant volatility σ i , we can obtain the dynamic process of an underlying asset with constant volatility as follows: The correlations between dW i (t) are defined by Equation (3).Jiang [20] carefully studied the explicit expression for a European multi-asset option value with constant volatility.Denote V BS (t, S(t); σ, Γ) as the corresponding price at time t, where the volatility vector is However, an analytic solution exists only for some specific options [20], such as exchange options, outer performance options, spread options, two dimension chooser options, basket options with a geometric average price, and quanto options.We give the specific expression for some of these in the numerical tests.

Dimension and Variance Reduction
In this section, we apply the acceleration methods of Liang and Xu [63] to price European multi-asset options with stochastic volatilities.The idea is that a martingale CV can be combined with the CMC method to reduce the variance of an MC simulation.

CMC Method
CMC can be used to reduce the variance of an MC simulation.Willard [64] initially put forward the CMC simulation to price options with stochastic volatilities.His method is feasible for those options that have a closed-form solution under the constant volatility model.Drimus [65] used CMC to analyze the variance products under the log-Ornstein-Uhlenbeck (log-OU) model.Boyle et al. [66] also used the CMC approach in pricing a down-and-in call option with a discretely monitored barrier.Broadie and Kaya [67] applied the CMC to accelerate exact simulations of the stochastic volatility with affine jump diffusion processes.Yang et al. [68] employed the CMC to reduce the variance of MC simulations when calculating the prices and greeks of barrier options.Dingeç and Hörmann [44] and Sun and Xu [46] combined CMC simulations and other variance reduction techniques to price basket options.
When we consider computing the expectation and the variance decomposition formula [69]: which indicates that the variance of E[Y|X] is always smaller than that of Y.The so-called CMC method uses the conditional expectation of the random variable E[Y|X] instead of that of the original random variable Y, which can obviously reduce variance.The key is that we need to have a closed form of E[Y|X].Now, we intend to deduce the CMC pricing formula for the pricing problem of Equation ( 8).The most important thing is to obtain the conditional expectation of the discounted payoff e −rT h(S(T)) on other random variables or stochastic information.First, a Cholesky decomposition of Brownian noise dW i (t) is conducted according to Equation (4): where dZ i (t) and d Zi (t) are independent standard Brownian noises, which means that: If we denote the vectors dW (t) = (dW 11) can be rewritten to: where n , and diag(v) is a diagonal matrix, with diagonal values from the vector v.According to Equations ( 12) and ( 13), it is obvious that: Notice that Equation (3) implies: We can solve for the covariance of d Z(t) by Equations ( 14) and ( 15) as: The entries of matrix Γ are To ensure the matrix Γ is well-defined, the condition of correlation coefficients Then, applying the Itô formula to ln(S i (t)), with the help of Equations ( 1) and (11), results in: ) Integrating both sides of the above equation from t to T results in: where Note that ξ i (t, T) is actually an exponential martingale with expectation E[ξ i (t, T)] = 1.Let σi (t, T) denote the average volatility of underlying asset S i on the interval [t, T], which is given by: It is observed that, given the information of stochastic processes {Z 1 (t), Z 2 (t), • • • , Z n (t)}, the quantities ξ i (t, T), and σi (t, T)(i = 1, 2, • • • , n) are totally determined.The pricing problem of Equation ( 8) then becomes the expectation of random variables { Z1 (t), Z2 (t), • • • , Zn (t)}.Assume there exists a Black-Scholes formula with constant volatilities.Then, calculating expectations with { Z1 (t), Z2 (t), • • • , Zn (t)} gives us the CMC pricing formula of a European multi-asset option with stochastic volatility, as follows: where • is the dot product of two vectors, ξ(t, T) = (ξ Compared with the MC formula in Equation ( 8), we now only need to simulate the n random variables {Z(t)} instead of the 2n random variables {W (t), and Z(t)}.Thus, the dimension and variance are reduced by the properties of the CMC.

Martingale Control Variate (CV)
To further reduce the simulation variance of the variable V BS t, S(t) • ξ(t, T); q • σ(t, T), Γ in Equation ( 22), a general martingale CV is proposed to combine with the CMC simulation.Some brief introductions about the CV method are given at first (for more details and references, please refer to Glasserman [39]).
When the CV method is used to compute the expectation E[Y], the CV estimator is: where X is called a CV and E[X] is the closed-form expectation of X.The constant b can be selected as b * = cov(X, Y)/var(X) to minimize the variance of the CV estimator with an optimal variance reduction ratio ).The efficiency of the CV method can be measured by the variance reduction ratio R 2 or the standard error reduction ratio R. The success of the CV depends on high correlations with the naive variable Y and the availability of its expectation E[X].Thus, the CV is always elegantly designed to a specific problem.In this paper, a martingale whose expectation equals to zero is used as a CV, which avoids any extra effort needed to obtain its expectation.
To construct an efficient CV in the CMC framework, a martingale representation theorem for the stochastic volatility pricing model of Equations ( 1) and ( 2) is proved in the following theorem.
Theorem 1 (Martingale Representation Theorem).If the underlying assets satisfy Equations ( 1) and ( 2), the European multi-asset option price at the initial time V(0, S(0), Y(0)) with payoff h(S(T)) can be expressed as: which can be rewritten in vector form as: where Proof.Applying Itô's formula to e −rt V(t, S(t), Y(t)) yields: where Furthermore, LV = 0 because of the Feynman-Kac formula [69], and thus: Now, by integrating both sides of the above equation on the interval [0, T], and noticing that V(T, S(T), Y(T)) = h(S(T)), we obtain the conclusion of the martingale representation theorem.
The martingale representation theorem gives us inspiration to construct efficient CVs.For simplicity, denote: The martingale expression in Equation (23) indicates that the variance of e −rT h(S(T)) in the MC simulation is totally determined by the martingale X.Thus, the martingale X plays the role of a perfect CV for an MC simulation.Fouque and Han [18] actually gave a similar representation in their work, and used the martingales as a CV to price single-asset options under a specific multi-factor stochastic volatility model.This can be understood in financial terminology.The martingale CV corresponds to a continuous delta hedge strategy taken by a trader who sells the option.The integrands of the martingale would, in theory, be the perfect delta hedges.Even though perfect replication by delta hedging under stochastic volatility models is impossible, the variance of replication error is directly related to the variance induced by the martingale CV method.
For the CMC pricing framework, taking the conditional expectation of both sides of Equation ( 23) based on the information {Z(t), 0 ≤ t ≤ T}, results in: where We can determine the expression for X by first substituting the Cholesky decomposition, Equation (11), into the expression of X in Equation (25).Then, we compute the expectations about where Equation (26) shows that the variance of V BS t, S(t) • ξ(t, T); q • σ(t, T), Γ is totally determined by the zero martingale X, since V(0, S 0 , Y 0 ) is a constant.This indicates that X, theoretically, is a perfect CV for V BS t, S(t) • ξ(t, T); q • σ(t, T), Γ in CMC simulations.It is a pity that we have no explicit expression of this perfect zero martingale, since there is no exact expression for V(t, s, y).A possible solution is that we approximate the option price V(t, s, y) with a Black-Scholes option price along with some carefully selected volatility parameters.In the following, we show our approach.
Given the information {Z(s), 0 ≤ s ≤ t}, the conditional expectation of S(t) can be computed by Equation (19) as: We intend to use V BS (t, S(t); θ, Γ) as an approximation of V(t, S(t), Y(t)), where θ is a constant vector whose values should be carefully selected.The partial derivatives can thus be approximated as: Now, given the information {Z(s), 0 ≤ s ≤ t}, by using the approximated derivatives and substituting S(t) with Ŝ(t) in Equations ( 29) and ( 30), F 1 (t) can be approximately expressed as: Notice that ∇ Y V BS (t, Ŝ(t); θ, Γ) = 0, so F2 (t) ≈ 0.
In the end, we obtain our martingale CV: The value of the constant volatility vector θ should be determined if we want to use the martingale X as a CV.Fouque and Han [18] illustrated a method for pricing a single-asset option with multi-factor volatility.They picked the long-term mean of the volatility as the volatility parameter in their specific multi-factor model.However, they did not offer a solution for non-mean-reverting stochastic volatility models, such as the Hull-White model.
In this paper, we set parameter θ = f (E[Y(t)]).The idea is that, on the interval [t, T], the stochastic variance Y(t) is approximated by the expectations of their initial state E[Y(t)].This results in a corresponding approximated stochastic volatility f (E[Y(t)]).We hope that the dynamic behavior of the approximated process with such parameters is similar to the original process.

Remark 1.
It is difficult to use the CMC pricing formula of Equation ( 22) if the analytic solution of a European multi-asset option price under constant volatility does not exist.However, we can still construct a martingale CV for an MC simulation in those cases.
According to the martingale representation theorem (Theorem 1), the variance of e −rT h(S(T)) in an MC simulation is totally determined by the zero martingale X (see Equation ( 23)).We can select a value V approx (t, S(t); θ, Γ) with a constant volatility parameter θ as the approximation of option price V(t, S(t), Y(t)) under a stochastic volatility model.Following the idea in the CMC framework, we select θ = f (E[Y(t)]).Thus, the partial derivatives can be approximated as: Furthermore, the martingale CV is: Taking the arithmetic average basket option with stochastic volatilities as an example, we can use the geometric average basket option with constant volatilities as an approximation and then construct the corresponding CV.It is expected that, for a more accurate approximated price V approx (t, S(t); θ, Γ), a larger variance reduction ratio can be obtained by the corresponding martingale CV.

Numerical Tests
In this section, we present some numerical tests designed for the typical multi-asset optionsincluding the exchange options, basket options and quanto options-to emphasize the efficiency of our method.

Exchange Options
The exchange option, which was first studied by Margrabe [70], empowers its holder with the right to exercise it by comparing the difference between the prices or the rates of return of two underlying assets.Its payoff is: If the underlying assets evolute with constant volatilities σ 1 and σ 2 , the exchange option has a pricing formula at time t, as shown by Margrabe [70] and Jiang [20] as follows: where and 2 dt is the cumulative distribution function of a standard normal variable.It is easy to derive the derivatives: We assumed the stochastic volatilities obey the Heston model, for i = 1, 2: The parameters should satisfy the Feller condition [71] to guarantee the positiveness of variance, i.e., 2a 1 θ 1 > σ 2  1 and 2a 2 θ 2 > σ 2 2 .We used a truncated Euler discrete scheme [39,72,73] with equal time intervals to simulate the Heston process in our tests.
At first, we wanted to examine the acceleration effect of a CMC simulation compared with a traditional MC simulation.We fixed the parameters S 1 (0) = S 2 (0) = K = 30, r = 0.05, T = 1, 015, and θ 2 = 0.05.Additionally, we took the number of time steps N = 100, and the number of simulations m = 100,000 in all numerical simulations.Note that < 1 should be satisfied from Equation (17).
Taking ρ 12 = 0 for simplicity, the numerical results with different correlation coefficients are recorded in Tables 2 and 3. Notes: ρ 1 is the correlation coefficient between the first underlying asset and its volatility; ρ 2 is the correlation coefficient between the second underlying asset and its volatility; V MC is the estimated price for MC simulation; and V CMC is the estimated price for CMC simulation. is the ratio of standard errors.
Table 2 records the estimated option values calculated by the MC and CMC simulations, which are denoted as V MC and V CMC , respectively.The upper part of Table 3 records the standard errors of an MC simulation, denoted as Std MC .The standard errors are almost the same for various correlation coefficients, and increase slightly with correlation ρ 2 while decreasing with ρ 1 .The exchange option can be seen as a call option on asset S 2 for fixed S 1 ; a higher correlation ρ 2 implies a larger variation in the price of asset S 2 , thus resulting in a larger value of the option price and a larger simulation variance.Similar analysis can be conducted with respect to correlation ρ 1 by regrading the exchange option as a put option on asset S 1 for fixed S 2 .
The middle part of Table 3 records the standard errors of a CMC simulation, denoted as Std CMC .Obviously, the standard CMC errors are always smaller than MC.It is interesting that a standard CMC error rapidly declines as correlation coefficient ρ 1 or ρ 2 tends to zero.Thus, the ratio of the standard errors of a CMC simulation to an MC simulation reduces.We denote this ratio as R = Std MC /Std CMC and present its values at the bottom of Table 3. R becomes larger when the correlation coefficient is getting closer to the original point, and decays rapidly in the opposite direction.For example, for ρ 1 = ρ 2 = 0, 0.25, 0.5, and 0.75, the reduction ratios of the standard error are 13.0340, 3.8805, 2.1170, and 1.4001, respectively.This can be explained by Equation (11); the CMC simulation removes the randomness that is independent from the stochastic variances Y 1 , Y 2 , and its quantity is proportional to In other words, a larger variance reduction ratio is promised when the absolute value of ρ 1 or ρ 2 is smaller.This property indicates that a CMC simulation is more competitive when the correlation between the underlying asset and stochastic volatility is weak.
We also investigated the computational costs of the MC and CMC methods.The computational platform for this paper was an Intel i5-6200U CPU, 2.30 GHz, 8 GB memory, and the software environment was Matlab R2018a for Windows 10.It took 50.88 s to calculate all of the data in the upper part of Table 3 and 25.85 s for the middle part, which means that the time cost of a CMC simulation is almost half that of an MC simulation.This is because the MC method needs to simulate four random variables, {W 1 (t), W 2 (t), Z 1 (t), and Z 2 (t)}, while the CMC method only needs to simulate two random variables, {Z 1 (t), and Z 2 (t)}.
Taking the variance reduction ratio into consideration, the speed up ratio of a CMC simulation to an MC simulation is defined as Thus, when correlation ρ 1 = ρ 2 = 0, the speed up ratio of the CMC is 13.0340 2 • 50. 88 25.85 = 334.38.Even for the case of a larger correlation ρ 1 = ρ 2 = 0.75, the speed up ratio of the CMC is 1.4001 2 • 50.88 25.85 = 3.86, which improves the efficiency of the MC simulation by roughly 75%.In summary, the CMC simulation enjoys the advantages of saving time and having a great variance reduction ratio, especially when the correlation coefficients are small.
We next tested the efficiency of our martingale CV method.As a contrast, we constructed another CV for the stochastic model, as suggested by Ma and Xu [74].Consider dummy assets whose prices Si (t), (i = 1, 2, • • • , n) satisfy the following stochastic differential equations: where σi (t) is a determined function.The covariance of dW i (t) is given by Equation (3).It can be computed by matching the first two moments of the underlying asset prices as σ2 In the case of a Heston stochastic volatility model: We used the payoff h excha ( S1 (T), S2 (T)) as a CV to the MC method, and we called this CV method a function CV method.The corresponding exchange option price can be computed using Equation (35) by replacing σ 2 1 + σ 2 2 − 2ρ 12 σ 1 σ 2 in Equations ( 36) and (37) with the average volatility on the interval [0, T] given by: 1 We changed the values of the correlation coefficients and kept the other parameters fixed as before.
Remember that The detailed results are shown in Table 4. Notes: ρ 12 is the correlation coefficient between the first underlying asset and the second underlying asset; ρ 1 is the correlation coefficient between the first underlying asset and its volatility; ρ 2 is the correlation coefficient between the second underlying asset and its volatility; Std MC , Std Mar , Std Fun are the standard errors from the MC method, the martingale CV method and the function CV method respectively; R 1 = Std MC /Std Mar ; and R 2 = Std MC /Std Fun .
In Table 4, Std MC , Std Mar , and Std Fun are the standard errors from the MC simulation, the martingale CV method, and the function CV method, respectively.R 1 = Std MC /Std Mar is the standard error reduction ratio of the martingale CV method compared to the MC simulation, and R 2 = Std MC /Std Fun is the the standard error reduction ratio of the function CV method compared to the MC simulation.
It is obvious that the standard error reduction ratio of the CMC is much larger than that of the function CV method, the former falling in 9-20 while the latter being about 3 or 4. Table 4 also shows that, for a fixed ρ 1 = ρ 2 , the standard errors of the MC simulation, martingale CV method, and function CV method decrease with the correlation value of ρ 12 .For a fixed ρ 12 , the standard errors of the MC simulation and function CV method increase with the value of ρ 1 = ρ 2 while the martingale CV method decreases with the absolute value of ρ 1 = ρ 2 , which is mainly caused by the properties of the CMC.Thus, the standard error reduction ratio of the martingale CV method also decreases with the absolute value of ρ 1 = ρ 2 .
The computing times for all values of the MC, the martingale CV, and the function CV methods are 22.33, 22.26, and 25.50 s, respectively.The time costs of the MC method and the martingale CV method are almost the same, while the function CV method is slightly slower.Thus, the martingale CV method proposed in our paper is superior to the function CV method, when considering the variance reduction ratio and the time cost.
Fixing the parameters ρ 12 = 0, and ρ 1 = ρ 2 = 0.5, we next examined the effects of the volatility parameters for the stochastic volatility.In the Heston stochastic volatility model, the Feller condition should be satisfied [71]; thus,  As shown in Table 5, the standard errors of the three simulation methods all increase with increasing volatilities of stochastic volatilities.Standard error reduction ratios also decline with the volatility of the stochastic volatilities.However, our martingale CV method is much more efficient than the function CV method, especially in the case of large volatility.

Basket Options
The payoff of the basket option at maturity depends on the average price of the underlying assets.Since the basket option with arithmetic average price does not have a closed-form price, even with constant volatility, we considered the geometric average basket option whose payoff at time T is: where n is the number of underlying assets, α i ≥ 0 are the weights of each underlying asset with ∑ n i=1 α i = 1, and K is the strike price.
The geometric average basket option has a closed-form solution if the underlying assets have constant volatilities as σ 1 , • • • , σ n .Denote: The geometric average basket option price at time t is given by Jiang [20]: where Thus, the derivatives are: For a basket option with n underlying assets, we still used the Heston stochastic volatility model and function CV method as a comparison.The expectation of the corresponding CV can be calculated using Equation (39) by substituting σ i σ j with 1 We fixed the parameters r = 0.05, T = 1, K = 30, S i (0) = 30, δ i = 0, a i = 2, and σ i = 0.2 (i = 1, 2, • • • , n).We allocated equal weights for the underlying assets, which means that α i = 1/n.For the initial value of the stochastic volatility, we took a linear interpolation between 0.1 2 and 0.3 2 for the n assets.In other words, the initial variance vector was Y(0) = (0.1 2 , 0.3 2 ) for n = 2 and Y(0) = (0.1 2 , 0.15 2 , 0.2 2 , 0.25 2 , 0.3 2 ) for n = 5.We took the long-term mean of stochastic variance as θ i = ( Y i (0) + 0.05) 2 , which was θ = (0.15 2 , 0.35 2 ) for n = 2, for example.For the correlations between Brownian noises, we took ρ ij = ρ 0 (i = j) for simplicity.To guarantee the positive definiteness of the matrix Γ = (ρ ij ), the parameter ρ 0 should satisfy −1/(n − 1) < ρ 0 < 1.In addition, < 1 is needed for the proper definition of Γ.Thus, we set ρ ij = ρ 0 = 0(i = j) at first.We fixed the number of time steps to N = 100 and the number of simulations to m = 100,000.We tested the acceleration effects of the CVs for different correlation coefficients ρ i .Numerical results are shown in Table 6.Notes: n is number of underlying assets; and ρ i is correlation coefficient between the ith underlying asset and its volatility.
Table 6 again shows that Std Mar , the standard error of the martingale CV method, decreases as the correlation coefficients ρ i tends to zero, resulting in a greater standard error reduction ratio R 1 in those cases.For example, R 1 goes from 30 to 9 when |ρ i | goes from 0 to 0.75.On the other hand, the simulation error Std Fun and, thus, the corresponding reduction ratio R 2 of the function CV method are not sensitive to the correlation coefficient.The reduction ratio is around 5 in all cases.Considering the number of underlying assets, the reduction ratio of the martingale CV slightly decreases as the number of assets n becomes larger, except for the ρ i = 0 cases.For example, the reduction ratios of the martingale CV method are 9.2740, 7.6660, and 6.3077 for n = 2, 5, and 10, respectively, and for ρ i = −0.75.On the other hand, the ratios are 30.5615,110.3248, and 291.4214 for the ρ i = 0 case.As a contrast, the performance of the function CV method is more stable with different n.It is obvious that our martingale CV is much more efficient than the function CV method.
Next, we fixed ρ i = 0.5(i = 1, 2, • • • , n) and changed the value of σ i , the volatility of the stochastic volatility.For convenience, we took an equal σ i for every underlying asset.The results are recorded in Table 7.As shown in Table 7, the standard errors of the three simulation methods increase with the volatilities of the stochastic volatility at fixed n, and decrease with the number of underlying assets for a fixed σ i .The standard error reduction ratios of the two CV methods decrease with increasing volatility of the stochastic volatility and increasing number of assets.However, the martingale CV method is more robust for different volatilities compared to the function CV method.For example, for the case of n = 2, the standard error reduction ratio of the martingale CV method decreases from 13.7479 to 12.2580 when σ i increases from 0.1 to 0.4, while that of the function CV method sharply decreases from 15.6830 to 4.1894.The results suggest that our martingale CV method is especially efficient in high volatility cases, while the function CV method has some advantages in a low volatility environment.

Quanto Options with Real Data
The quanto option is a contract written when someone invests money in foreign securities.Usually, its risk depends on the volatility of the securities' prices and the change of the foreign currency rate.Its main purpose is to provide exposure to a foreign asset without taking the corresponding exchange rate risk.We applied our method to price a quanto option.Park et al. [62] used a power series expansion method to obtain an analytic approximation value for the quanto option price under the Hull-White stochastic volatility model.
In Table 9, ρ 12 stands for the correlation between the S&P500 and FX rate.V Appro is the approximated value obtained by the series expansion method in Park et al. [62].V MC , V Mar , and V Fun are the estimated values of the MC simulation, the martingale CV method, and the function CV method, respectively.Std MC , Std Mar , and Std Fun are the standard errors of the MC simulation, the martingale CV method, and the function CV method, respectively.R 1 = Std MC /Std Mar is the standard error reduction ratio of the martingale CV method compared to the MC simulation, and R 2 = Std MC /Std Fun is the the standard error reduction ratio of the function CV method compared to the MC simulation.For the function CV method, the expectation of the corresponding CV can be calculated by using Equation ( 41) and substituting σ i σ j with 1 T T 0 E[σ i (t)]E[σ j (t)]dt, where E[Y i (t)] = Y i (0)e µ i t , i = 1, 2. It is obvious that our martingale CV method has a larger standard reduction ratio than the function CV method.This, again, shows the efficiency and robustness of our method.

Conclusions
In the context of European multi-asset options with stochastic volatilities, we propose a dimension and variance reduction Monte Carlo method.A conditional Monte Carlo pricing formula is deduced, and then the martingale representation theorem is proved.A martingale control variate is combined with the conditional Monte Carlo simulation.
Numerical tests on typical multi-asset options-including exchange options, basket options, and currency options-showed that this method yields considerable variance reduction, not only when compared to a traditional Monte Carlo simulation, but also with respect to the function control variate in Ma and Xu [74].
For future research, it would be interesting and challenging to extend the framework in this paper to price more options with stochastic volatilities, not only European options but also exotic options such as American options or barrier options.Furthermore, it would be interesting to study jump diffusion models with stochastic volatilities.Another important approach is to use this framework in empirical financial studies and risk management.After model parameters are calibrated with real market data, our method can be used to accurately and quickly value option prices which can be widely used in areas of economics and finance.We would also like to extend this method to other areas like risk management and civil engineering.

Table 1 .
Models of stochastic volatility.

Table 2 .
Exchange option: Estimated prices for MC and CMC with different correlation coefficients.

Table 3 .
Exchange option: Numerical results for MC and CMC with different correlation coefficients.1 is the correlation coefficient between the first underlying asset and its volatility; ρ 2 is the correlation coefficient between the second underlying asset and its volatility; Std MC , Std CMC are the standard errors of estimated prices from the MC and CMC simulations respectively; and R = Std MC

Table 4 .
Exchange option: Numerical results for CVs with different correlation coefficients.

Table 5 .
Exchange option: Numerical results for CVs with different volatilities of the stochastic volatilities.

Table 6 .
Geometric average basket option: Numerical results for CVs with different correlation coefficients.

Table 7 .
Geometric basket option: Numerical results for CVs with different volatilities of the stochastic volatility.
Notes: n is number of underlying assets; and σ i is the volatility of the ith stochastic volatility.

Table 9 .
Quanto option: Numerical results with different correlation coefficients.