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