How to do a Brownian Bridge with quasi-random numbers in the Heston model?

Question

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
    %variance.

    %Generate correlated random numbers

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

    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, :);
    end

    option = exp(-r*T)*max(exp(X(end, :))-K, 0);
    price = mean(option);
    error = std(option)/sqrt(m);
end

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$

Jianwei Zhu's example

Output