Numerical Simulation of the Heston Model under Stochastic Correlation

.


Introduction
The Heston Model Heston (1993) is one of the most widely used stochastic volatility models.It is an extension of the Black-Scholes (Black and Scholes 1973) model by taking into account stochastic volatility given by the Cox-Ingersoll-Ross (CIR) process.The attractiveness of the Heston model is its analytical tractability and the consideration of the correlation between the underlying asset price process and volatility process.Subsequently, a couple of papers on the numerically stable and efficient computation of European-style option prices were published, e.g., (Andersen 2008;Carr and Madan 1999;Kahl and Jäckel 2005;Lee 2004;Lewis 2001;Lipton 2002).
It has been pointed out, in many works (see, e.g., (Christoffersen et al. 2009;Grzelak and Oosterlee 2011)) that the Heston model is unable to provide enough skew in the implied volatility as market required, especially for a short maturity.Therefore, one tries to extend the Heston model.For example, one way is to extend the Heston by introducing a more realistic stochastic volatility process, which is the double Heston model (Christoffersen et al. 2009) or by introducing a stochastic interest rate, which is the Hybrid-Heston-Hull-White model (HHW) (Grzelak and Oosterlee 2011); another way is to adapt the Heston model by allowing time-dependent parameters (see Elices 2009;Mikhailov and Nögel 2003).Actually, the major factor that affects the implied volatility skew is the correlation.However, in the pure Heston model (Heston 1993), and also in most of the extended Heston models, only a constant correlation coefficient is used.It has been shown in Teng et al. (2014bTeng et al. ( , 2016b) ) that the calibration to market data can be already improved by allowing a local time-dependent correlation model.Therefore, in this work, we extend the Heston model by imposing a stochastic correlation model cf.Emmerich (2006); Teng et al. (2014aTeng et al. ( , 2016a) ) and discuss its simulation methods.Furthermore, using Monte Carlo simulation, we study how the implied volatility is affected by introducing stochastic correlations.
In the next section, we impose a stochastic correlation model to the Heston model and discuss in Section 3 a discretization for each path of the variance, correlation and the log price process including numerical analysis.Section 4 is devoted to firstly comparing different numerical algorithms, and secondly to recognizing the effect of imposing stochastic correlation on implied volatility.Finally, Section 5 concludes this work.

Stochastic Correlation in the Heston Model
We study the Heston model and its extension by incorporating a stochastic correlation process.Heston's stochastic volatility model under the risk-neutral measure reads dS t = rS t dt + √ ν t S t dW S t , (1) where S t denotes the spot price of the asset and ν t is the instantaneous variance, where µ ν is the long-term variance, κ ν is the speed at which it reverts to µ ν and σ ν is the volatility of the variance process.We note that the process ( 2) is strictly positive if the parameters obey the Feller condition 2κ ν µ ν > σ ν .The Brownian motions (BMs) W S and W ν are correlated with a constant ρ Sν by and under risk-neutral measure.By applying Itô's Lemma with x t = log(S t ) (log-transform), we obtain from (1) the log price process as where W x t is the same BM to W S t in (1).We suppose an appropriate stochastic correlation process of the form where ã(t, ρ t ) and b(t, ρ t ) are given functions depending on the chosen correlation process, Wρ t is a standard BM with respect to the physical measure.By including the market price of correlation risk, the correlation process (5) can be rewritten as which is under risk-neutral measure, where λ(S t , ν t , ρ t , t) represents the price of the correlation risk and could be assumed to be λρ t , with a constant λ.In what follows, we set ã(t, With the aim of imposing a stochastic correlation between the log price process dx t and the stochastic variance dν t , we extend the Heston model as with where ρ t is given by (8), and ρ xρ and ρ νρ are assumed to be two constant correlations.

Path Simulation
We now discuss how to simulate the paths for ( 7)-( 9) to compute the price of European options in the extended Heston model applying Monte Carlo simulation.We need to generate random paths of the triplet (ν t , ρ t , x t ) for all t ∈ {t i } N i=1 := T .To be more precise, for an arbitrary time increment ∆, we need to generate a random sample of (ν t+∆ , ρ t+∆ , x t+∆ ) for given (ν t , ρ t , x t ).Repeated application of the resulting one period scheme will generate a full path (ν t , ρ t , x t ) t∈T .

Discretization for the Variance Process ν t
To discretize the variance process ν t , we employ the quadratic-exponential (QE) scheme by Andersen (2008).The idea of the QE scheme is the approximation to the non-central chi-square distribution.Let νt denote a discrete-time approximation to ν t , for sufficiently large realized values of νt , Andersen (2008) suggested to approximate the non-central chi-square random variable by the power function νt+∆ = α(β where Z ν is a standard Gaussian random variable, α and β are certain constants that can be determined by moment-matching using the parameters in (7) and given in the following Proposition 1; the detailed calculations can be found in Andersen (2008).
Proposition 1.The mean and the variance of the variance process (7) read If we set ψ := s 2 m 2 and choose and then (11) has a mean equal to m and a variance equal to s 2 .Note that ψ ≤ 2.
However, the approximation (11) will not work well for small values of νt , cf.Andersen (2008).For small values of νt , Andersen (2008) suggested to use an approximated density for νt+∆ of the form: where δ is a Dirac delta function, and p ∈ [0, 1] and q > 0 are constants to be determined by moment-matching.We integrate ( 16) and obtain and, by inverting, we obtain Thus, the sampling scheme for small values of νt reads where U ν is a uniform random variable.
Proposition 2. Let m, s 2 and ψ be defined as in Proposition 1.For ψ ≥ 1, there exist two parameters p and q such that (19) has a mean equal to m and a variance equal to s 2 , which read We only need to select an arbitrary level ψ c ∈ [1, 2] and choose either (11) or ( 19) according to ψ ≤ ψ c or ψ > ψ c to do the sampling for the variance process (7).

Discretization for the Correlation Process ρ t
We know that several stochastic processes could be used for modelling stochastic correlation, e.g., the bounded Jacobi process (Ma 2009;Emmerich 2006), the sort of stochastic correlation process produced by transformation (Teng et al. 2014a(Teng et al. , 2016a)).Moreover, for nice analytical tractability, we might choose some other mean-reverting processes that exhibit simpler structure, e.g., the Ornstein-Uhlenbeck (OU) process (Uhlenbeck and Ornstein 1930).
with the exact solution where Z ρ is a standard Gaussian random variable.Thus, the functions a(t, ρ t ) and b(t, ρ t ) defined in (5) are known as κ ρ (µ ρ − ρ t ) and σ ρ , respectively.However, the major drawback of using the OU process for stochastic correlation is that the process is not bounded.This is to say the generated correlations can be out of the correlation interval [−1, 1], especially with a small value of κ ρ and a large value of σ ρ .Moreover, it has been indicated by Teng et al. in (Teng et al. 2016b) that and the condition for P( This does not necessarily means that σ ρ tends to zero.If one limits the mean value µ ρ to be in (−1, 1), from ( 24) and ( 25), one can conclude that the OU process is bounded in the interval with the condition In practice, such a positive constant C could be selected such that the condition σ ρ ≥ C is already sufficient to ensure that the generated correlations stay in (−1, 1), if the initial correlation ρ 0 and the long-term mean µ ρ are not close to −1 or 1.

Discretization for the Log Price Process x t
In this section, we discuss how to discretize the log price process (9).As indicated by Andersen (2008), a straight discretization of (9) may lead to the problem of "leaking correlation": suppose that we use directly a Euler scheme for simulating (9): We know that the true correlation between xt+∆ and νt+∆ is always close to ρ t given by ( 8).However, νt+∆ and Z ν in (11) have a strong nonlinear relationship, which will imply that the effective correlation between xt+∆ and νt+∆ will be closer to zero than ρ t for the cases where the probability P(β + Z v < 0) is nonzero.To tackle this problem of "leaking correlation", we reformulate the stochastic differential equation (SDE) system ( 7)-( 9) as follows: first, we see that the SDE system ( 7)-( 9) has a family of correlation matrices To simplify notation, we set ρ 1 := ρ νρ (ρ ρν ) and ρ 2 := ρ xρ (ρ ρx ).One can thus perform a Cholesky-decomposition C t = L t L t , where L t is a family of lower triangular matrices given by which can be used to reformulate the SDE system ( 7)-( 9) with respect to the independent BMs Wν t , Wρ t and Wx t as: Since our main aim is to impose a stochastic correlation between the asset process (31) and the stochastic variance process (29), to simplify the model, we assume ρ 1 = 0, and the latter SDE system thus becomes In the following, we will discuss two different ways to discretize (34).

The Euler and Milstein Scheme (EM Scheme)
The most simple way of discretizing ( 34) is to use the Euler or Milstein scheme.The discretization of (34) by applying the Euler scheme (Euler [1768(Euler [ ] 1913;;Maruyama 1955) where Z ρ, Z ν and Z x are independent standard Gaussian random variables.The discretization of (34) by applying the Milstein scheme Milstein (1974) will be the same to ( 35), since all the derivatives included in the coefficients of the double integral terms (with respect to BMs) by the Milstein scheme are equal to zero.Moreover, we remark that Ŝt = exp( xt ) with the discretized process xt in ( 35) is a martingale, and any types of stochastic correlation processes can be straightforwardly employed within this scheme.

The Hybrid Scheme (HB Scheme)
In this section, we introduce a new way to discretize (34) where several different approximation techniques will be used.We thus call it the hybrid scheme.We take the OU process as an example due to its analytical tractability and start with the integral form of (34) where ρ t follows an OU process.For the integral , where the integrand is not independent from the corresponding BM, we consider Itô's product in the following where it has been assumed that ν t and ρ t are independent from each other corresponding to ρ 1 = 0.This product implies that Now, we insert (38) into (36) For the integrals over the time in (39), we simply use the approximation and where γ 1 and γ 2 are given constants, e.g., we choose γ 1 = γ 2 = 1 2 for a central discretization.In all the Itô integrals in (39), the integrand is independent with the corresponding BM.Thus, they can be approximated by and respectively, where Z ρ and Z x are independent standard Gaussian random variables.With all the approximations mentioned above, we rearrange (39) as where A analysis of the convergence properties for ( 44) is difficult and complicated, since it may not have any high-order moments.We will consider the analysis of weak consistency.
Proposition 3. Assume that γ 1 + γ 2 in (44) approach 1 for ∆ → 0. Conditional on Ŝt , νt and ρt , we have for the HB scheme lim lim Proof.The statements ( 46) and ( 47 The first part in (45) can be proved in the same way.Moreover, we calculate In the same way, one can prove the second part in (48).
Obviously, Proposition (3) says the HB scheme is weakly consistent Kloeden and Platen (1999) whilst which will be satisfied with non-extreme parameter values.

HB Scheme with Martingale Correction (HBM Scheme)
We know that the price process S t will be a martingale; however, the price process S t = exp(x t ) in ( 44) is not a martingale.For this problem, on one side, we can reduce the size of ∆; on the other side, the "martingale correction" proposed by Andersen (2008) can be employed.
The scheme ( 44) is equivalent to Assuming that ρ t+∆ is known, by iterated conditional expectations, we calculate that Clearly, it is not easy to compute the part underlined in the latter equation.However, since both exponents of νt in that part are greater than one, we might thus ignore this part to obtain an approximated martingale correction.Therefore, we reformulate the latter equation as where We assume that M in (51) is finite, and then E Ŝt+∆ | Ŝt , ρt+∆ is also finite.In order to force E Ŝt+∆ | Ŝt , ρt+∆ will be a martingale, we require which implies K 0 = − ln M − N. Finally, we obtain the HBM scheme as Obviously, the challenge of using the HBM scheme is to compute K 0 , which is actually a random variable conditional on Ŝt and ρt+∆ .Because ν t and ρ t are independent, hence, we can directly adopt the recent approach by Andersen (2008, Proposition 9) where α, β, ψ, ψ c , p, q are defined in Propositions 1 and 2, and N and A are defined in ( 52) and ( 53).

Numerical Results
As mentioned before, we analyze the numerical results by pricing European call options.We denote the exact option price C and the numerically approximated price with Ĉ that can be computed using the expectation and approximated by a Monte Carlo method Thus, we define the error of a discretization scheme as which will be dependent on ∆.For all of the numerical experiments, we assume S = 100, r = 0 and three different levels of the strike K = [70, 100, 140].

A Comparison of the Numerical Methods EM, HB and HBM
In this section, we test the discretization schemes EM, HB and HBM described in Section 3.2 for the log price process x t .It is well known that the OU process is a mean-reverting process, i.e., whilst we initialize the stochastic correlation process so that it can rapidly reach its mean value µ ρ , the option price computed in the extended Heston model should be the same as the original Heston price with the constant correlation ρ = µ ρ .In this case, we can take the original Heston price obtained by computing the (semi-)analytical pricing formula in Heston (1993) as the benchmark.
To initialize the variance process, we take the parameters collection used in Andersen (2008) and given in Table 1.We see that all the parameters collection of Cases I, II, III are not under the Feller condition.Hence, for Case IV, we choose parameters collection that are under the Feller condition.In order to let the price be computed in the extended Heston model to coincide with the pure Heston price, as mentioned above, e.g., we choose κ ρ = 2, σ ρ = 10 −3 and set the value for µ ρ and ρ 0 to be same as the value of ρ in Table 1.Moreover, letting ρ 2 = 0, γ 1 = γ 2 in (40) be 0.5, we use M = 10 6 for the Monte Carlo method and report the errors for different time steps and for cases in Tables 2-5 by varying the value of the time step ∆ from one year to 1/32 year.We consider first the results for Case I in Table 2 and find that both the discretization schemes HB and HBM have an advantage over the EM scheme, and the advantage is considerable for the out-of-money options with K = 140.By comparison to the HB scheme, one can realize that it is beneficial by adding a martingale correction for computing the in-the-money and at-the-money options with a simulation step of ∆ = 1 or ∆ = 1 2 .Since the results for Cases II and III are qualitatively similar to those of Case I, one can reach the same conclusion that both the HB and HBM schemes outperform the EM scheme.It is worth to noting, for the less challenging Case III, that the results for Case III by using the EM scheme are better than those of Cases I and II.Now, we turn to the Case IV where the parameters obey the Feller condition.The performance of the HB scheme in this case is a bit poor, especially, with a large time step ∆.Actually, we have also tested the pure QE scheme in Andersen (2008) for this case and obtained the same output.By contrast, the performance of the EM scheme for this case is surprisingly good.Fortunately, for this case, the martingale correction has brought a huge benefit so that the HBM scheme still performs better.
About comparing the numerical efficiency, we check the average computation times of the HB and HBM schemes relative to the EM scheme for all runs in Tables 2-5, which are 0.61 and 0.79, respectively.Obviously, both HB and HBM schemes are more efficient than the EM scheme because, compared to the EM scheme, we have one random variable less to simulate with the HB and HBM schemes.

The Effect of Imposing Stochastic Correlation on Implied Volatility
In this section, we analyze the effect of imposing stochastic correlation on the implied volatilities.
To do this, we show the role of using stochastic correlation process in implied volatility, namely to see how the values of parameters of the correlation process will drive the implied volatilities.We display in Figure 1 the changes of implied volatilities by varying each parameter of stochastic correlation process.For this experiment, we prefer to use the HBM scheme.We consider a Call-option with S = 120, T = 0.5 and the strikes from 114 to 126 in increments of 1, r = 1%.For the variance process, we set ν 0 = 0.03, µ ν = 0.04, κ ν = 2.1, σ ν = 0.4, and, for the correlation process, we choose κ ρ = 3.5, σ ρ = 0.1, µ ρ = −0.6 (equal to the constant correlation) and set ρ 0 = −0.4except for the one that is varying.Finally, we set ρ 2 = 0.1 and use M = 10 6 for the Monte Carlo simulation.
From Figure 1, we realize that the parameters of correlation process can control the skewness or smiles of the implied volatilities.Compared to using a constant correlation parameter, including stochastic correlation provides more flexibility and can thus improve the calibration to the real market data.

A Comparison with the Effect of Stochastic Correlation
In order to be able to compare with the pure Heston price, in the numerical experiments above, we have initialized the stochastic correlation process such that it can rapidly reach its mean value.Now, we want to test our numerical schemes including the effect of stochastic correlations.Teng et al. (2016c) found the well approximated characteristic function of the Heston model extended by including stochastic correlations driven by the OU process in a closed-form, which can be used for analytical pricing purposes.Moreover, a comparison to the EM scheme has already been provided in Teng et al. (2016c); for more detailed information, we refer to Teng et al. (2016c).
Next, we test all of the numerical schemes by comparing them to the method in Teng et al. (2016c).For this, we take a five years Call with S 0 = 100 and report the price differences between using the proposed numerical schemes and the approach in Teng et al. (2016c) in Table 6.The used parameter values for dν t are not under the Feller condition.We thus obtain similar results to those of Case IV: the HB scheme does not perform well with a large time step ∆.From all the numerical experiments, we conclude that: the HB and HBM schemes outperform the EM scheme when the parameter values of dν are not subject to the Feller condition.By contrast, under the Feller condition, the EM and HBM schemes are more favorable since they outperform the HB scheme at a large time step ∆.

Conclusions
In this work, we extended the Heston model by imposing stochastic correlations driven by a SDE.We have introduced different numerical algorithms including numerical analysis and compared their merits.We showed which algorithms are more favorable for which model parameterizations.Thereof, the HB and HBM scheme arose by adopting the QE scheme in Andersen (2008).
A couple of numerical results are provided.It has been shown that the numerical schemes proposed in this paper can work so well for the extended Heston by including stochastic correlations as the QE scheme in Andersen (2008) for the pure Heston model.Moreover, we realized the benefit of incorporating stochastic correlations by investigating the effect of stochastic correlations on the implied volatility.Because of the increased number of model parameters through the correlation process, the extended Heston with stochastic correlations can provide more than enough skews or smiles in the implied volatility as the market requires than the pure Heston model.

Table 1 .
Parameters collection for dν t .

Table 2 .
A comparison of the relative errors in Case I using different schemes; numbers in parentheses are standard deviations.

Table 3 .
A comparison of the relative errors in Case II using different schemes; numbers in parentheses are standard deviations.

Table 4 .
A comparison of the relative errors in Case III using different schemes; numbers in parentheses are standard deviations.

Table 5 .
A comparison of the relative errors in Case IV using different schemes; numbers in parentheses are standard deviations.