I'm required to use the Euler Monte Carlo method to compute the option price under Heston model settings.

I know from some paper that the convergence is volatile for the Heston model with a plain Monte Carlo and Euler scheme, so I use the Sobol sequence to generate quasi-random numbers. However, now the problem comes to that the generated series is correlated at different time steps.

In detail, it's OK when I use norminv to transform the quasi-random numbers to standard normal distribution. However, after I use the correlation $\rho$ to make correlated series at the same time step, it turns out that the series in different time steps are also correlated. If I ignore this problem, the result of my function would be significantly larger than the analytical solution given by https://quanto.inria.fr/premia/koPremia.

Should I use a Brownian Bridge to deal with this problem? If so, how do I do it in MATLAB? Or how could I eliminate the series correlation in quasi-random numbers?

Here is my MATLAB code for the plain Monte Carlo method. The problem is that it doesn't converge at all.

function [price,error] = heston_mc(S0, K, V0, rho, kappa, theta, eta, r, T, n, m, rseed)

    %V0 is the initial variance, and kappa is the mean reversion speed of
    %variance, theta the mean reversion level, neta the volatility of

    %Generate correlated random numbers

    if ~exist('rseed','var')
        rseed = 42;

    seed = RandStream.create('mcg16807','seed',rseed);
    z = randn(seed, n-1, 2*m);
    zX = z(:, 1:m);
    zV = rho*zX + sqrt(1-rho^2)*z(:, m+1:2*m);
    delta = T/n;
    V = [repmat(V0,1,m); zeros(n-1,m)];
    X = [repmat(log(S0), 1, m); zeros(n-1, m)];

    for i = 2:n
        X(i, :) = X(i-1, :) + (r-max(V(i-1, :), 0)/2)*delta + sqrt(delta)*sqrt(max(V(i-1, :), 0)).*zX(i-1, :);
        V(i, :) = V(i-1, :) + kappa*(theta-max(V(i-1, :),0))*delta + eta*sqrt(delta)*sqrt(max(V(i-1, :), 0)).*zV(i-1, :);

    option = exp(-r*T)*max(exp(X(end, :))-K, 0);
    price = mean(option);
    error = std(option)/sqrt(m);
  • 3
    $\begingroup$ quasi randoms won't help with the problem that the Euler scheme is bad for Heston and requires an awful lot of time steps to get good behaviour. There have been a lot of papers on better schemes now. See ssrn.com/abstract=1617187 and references therein. $\endgroup$
    – Mark Joshi
    Commented Jul 5, 2015 at 22:27

1 Answer 1


You don't need to use the Sobol sequence to generate quasi-random numbers in MATLAB. We know the Heston model is represented by the bi-variate system of stochastic differential equations (SDE):

\begin{align} & d{{S}_{t}}=rS_tdt+{\sqrt\upsilon_t} d{{W}_{1}}(t) \\ & d{{\upsilon}_{t}}=\kappa(\theta-\upsilon_t) dt+\sigma{\sqrt\upsilon_t}d{{W}_{2}}(t), \\ \end{align} where \begin{align} dW_1(t)dW_2(t)=\rho dt. \end{align}

All of the simulation schemes, like the Euler-scheme for the Heston model, contain the same basic steps. First, two independent standard normal random variables are generated, and then made dependent by applying a Cholesky decomposition. These are multiplied by $\sqrt{\Delta t}$ to make them proxy Brownian motion increments. Second, we obtain the updated $v_{t+\Delta t}$. Third, we obtain the updated value $S_{t+\Delta t}$ or $x_t=\ln S_{t+\Delta t}$.

  1. Initialize $S_0$ to the spot price (or $x_0$ to the log spot price), and initialize $v_0$ to the current variance parameter.
  2. Generate two independent normal random variables $Z_1$ and $Z_2$, and define $Z_v= Z_1$ and $Z_S=\rho Z_v+\sqrt{1-\rho^2}Z_2$. Proxy the Brownian motion by $dW_1(t)=\sqrt{\Delta t}Z_S$ and $dW_2(t)=\sqrt{\Delta t}Z_v$.
  3. Obtain the updated value $v_{t+\Delta t}$.
  4. Given $v_{t+\Delta t}$, obtain the updated value $S_{t+\Delta t}$ (or $x_{t+\Delta t}$) and return to Step 2.

Implementation in the MATLAB

Please insert $\kappa, \theta, \sigma, v_0, S_0, r, q, \rho, T, M, and N$

Enter image description here

Jianwei Zhu's example

Enter image description here


Enter image description here

  • 2
    $\begingroup$ @ zackhui your euler sheme differs from the one of Behrouz, since your V is allowed to have negative startvalues in each step. maybe that's the problem. $\endgroup$
    – Phun
    Commented Jul 5, 2015 at 16:24
  • 2
    $\begingroup$ @Phun probably I'll try to change to a fully truncated method, thanks for pointing out. $\endgroup$
    – zackhui
    Commented Jul 5, 2015 at 17:50
  • 2
    $\begingroup$ @ zackhui you can choose two mthtods :$v_{t+\Delta t}\rightarrow |v_{t+\Delta t}|$ or $v_{t+\Delta t}\rightarrow max\{v_{t+\Delta t},0\}$ $\endgroup$
    – user16651
    Commented Jul 5, 2015 at 17:59
  • 2
    $\begingroup$ @BehrouzMaleki exactly and I choose the second one. $\endgroup$
    – zackhui
    Commented Jul 5, 2015 at 18:03
  • 2
    $\begingroup$ Hi both, I think there is a confusion between pseudo random numbers (e.g. generated using Mersenne-Twister algorithm which is the default used by MATLAB AFAIR) and quasi random numbers (low discrepancy pseudo random sequences e.g. Sobol). But this does not change anything to what you both sad. Just saying :) $\endgroup$
    – Quantuple
    Commented May 26, 2016 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.