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