Simulation of Heston process

I am currently working on implementing Heston model in matlab for option pricing (in this case I am trying to price a European call) and I wanted to compare the results I obtain from using the exact formula and the Monte-Carlo simulation using the Milstein discretization.

In order to provide you more details about what I am doing, here is the equation I use for each code :

The process :

$dS_t = \mu_tS_tdt + v_t^{1/2}S_tdZ_1$

$dv_t = -\lambda(v_t - \theta)dt + \eta(v_t)^{1/2}dZ_2$ with $<dZ_1 dZ_2> = \rho dt$

Exact formula :

In this part, I am using the following formula (from Gatheral). $C(x, v, \tau) = K[e^xP_1(x, v, \tau) - P_0(x, v, \tau)]$

$\alpha = -\frac{u^2}{2} - \frac{iu}{2} + iju$

$\beta = \lambda - \rho \eta j - \rho \eta i u$

$\gamma = \frac{\eta^2}{2}$

$r_\pm = \frac{\beta \pm \sqrt{\beta^2 - 4\alpha \gamma}}{2\gamma} := \frac{\gamma \pm d}{\eta^2}$

$g := \frac{r_-}{r_+}$

$D(u, \tau) = r_- \frac{1-e^{-d\tau}}{1-ge^{-d\tau}}$

$C(u, \tau) = \lambda[r_- \tau- \frac{2}{\eta^2}log(\frac{1-ge^{-d\tau}}{1-g})]$

$P_j(x, v, \tau) = \frac{1}{2} + \frac{1}{\pi}\int_0^{\infty}{Re(\frac{exp(C_j(u, \tau)\theta + D_j(u, \tau)v+iux)}{iu}du)}$

Monte-Carlo scheme

$v_{i+1} = (\sqrt{v_i} + \frac{\eta}{2} \sqrt{\Delta t}Z)^2 - \lambda (v_i - \theta)\Delta t - \frac{\eta^2}{4} \Delta t$

$x_{i+1} = x_i - \frac{v_i}{2}\Delta t + \sqrt{v_i \Delta t}W$ with $x_i := log(S_i/S_0)$

The problem I have when implementing this two methods is that the result I get is very different from each other (by factor 10). I compared the result of the first method with some I could find on internet and it was coherent with what I get. I also compare the result of the Monte-Carlo simulation with parameters with bring back the Black-Scholes formula and the result is quite good.

The problem arise when I try the same parameters on the 2 methods and the output is no even comparable from each other.

I provide you also with the code I use for implementing those methods.

Close form method :

function [ res ] = HestonPrice( t, S, V, K, T, r, lambda, vMean, eta, rho)
%UNTITLED Summary of this function goes here
%   This function compute the price of a call option using the close form
%   of the Heston model
%   Inputs :
%   t : date at which one we want to compute the price
%   S : price of the asset at time t
%   V : volatility at time t
%   K : Strike of the option
%   T : Expiration date of the option
%   r : Risk free rate
%   lambda : speed of the mean reversion
%   vMean : long-term mean for the variance (not the volatility !!!)
%   eta : volatility of the volatility equation
%   rho : correlation between the two brownian processes
%
%   Output :
%   res : Price of the call option
%
%   Example :
%   HestonPrice(0, 1, 0.16, 2, 10, 0, 1, 0.16, 2, -0.8) -> 0.0495

FtT = S * exp(r*(T-t));
x = log(FtT / K);
tau = T-t;

P0 = P(0, V, x, tau, rho, eta, lambda, vMean);
P1 = P(1, V, x, tau, rho, eta, lambda, vMean);

res = K * (exp(x) * P1 - P0) * exp(-r*(T-t));

end

function [ res ] = D(j, u, tau, rho, eta, lambda)

alpha = - u .* u / 2 - 1i .* u / 2 + 1i .* j .* u;
beta = lambda - rho .* eta .* j - rho .* eta .* 1i .* u;
gamma = eta .* eta / 2;
d = sqrt(beta .* beta - 4 .* alpha .* gamma);

rPlus = (beta + d)./(2 * gamma);
rMinus = (beta - d)./(2 * gamma);
g = rMinus ./ rPlus;

res = rMinus .* (1 - exp(-d .* tau)) ./ (1 - g .* exp(-d .* tau));

end

function [ res ] = C(j, u, tau, rho, eta, lambda)

alpha = - u .* u / 2 - 1i * u / 2 + 1i .* j .* u;
beta = lambda - rho * eta * j - rho .* eta .* 1i .* u;
gamma = eta .* eta / 2;
d = sqrt(beta .* beta - 4 .* alpha .* gamma);

rPlus = (beta + d)./(2 .* gamma);
rMinus = (beta - d)./(2 .* gamma);
g = rMinus ./ rPlus;

res = lambda .* (rMinus .* tau - 2 ./ (eta .* eta) .* log((1 - g .* exp(-d .* tau))./(1-g)));

end

function [ res ] = funToIntegrate(j, v, x, u, tau, rho, eta, lambda, vMean)
CCompute = C(j, u, tau, rho, eta, lambda);
DCompute = D(j, u, tau, rho, eta, lambda);
res = real(exp(CCompute .* vMean + DCompute .* v + 1i .* u .* x) ./ (1i .* u));
end

function [ res ] = P(j, v, x, tau, rho, eta, lambda, vMean)
fun = @(u) funToIntegrate(j, v, x, u, tau, rho, eta, lambda, vMean);
resIntegration = integral(fun, 0, Inf);

res = 1/2 + 1/pi * resIntegration;
end


Monte Carlo Method :

function [ res ] = HestonPrice_MC( t, S, V, K, T, r, lambda, vMean, eta, rho, n, samples )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
%   HestonPrice_MC(0, 1, 0.16, 2, 10, 0, 1, 0.16, 2, -0.8, 1000, 10000)

temp = zeros(samples, 1);

parfor i=1:samples
temp(i) = finalPriceX(t, r, V, T, lambda, vMean, eta, rho, n);
end
temp = S * exp(temp);
temp = max(temp - K, 0);
res = mean(temp) * exp(-r * (T-t));
end

function [X, V] = nextStep(r, x, v, lambda, vMean, eta, rho, deltaT)
coovMat = [1 rho; rho 1];
rdVar = mvnrnd([0 0], coovMat);

V = (sqrt(v) + eta * 0.5 * sqrt(deltaT) * rdVar(1))^2 - lambda * (v - vMean) * deltaT - eta^2 * 0.25 * deltaT;
if V < 0
V = 0;
end
X = x + (r - v * 0.5) * deltaT + sqrt(v * deltaT) * rdVar(2);
end

function [price] = finalPriceX(t, r, V, T, lambda, vMean, eta, rho, n)
x = 0;
v = V;
deltaT = (T-t)/n;

for i=1:n
[x_temp, v_temp] = nextStep(r, x, v, lambda, vMean, eta, rho, deltaT);
x = x_temp;
v = v_temp;
end

price = x;
end


I use the following parameters for the simulation :

fprintf('Computation of the Heston price with the close form formula');
HestonPrice(0, 1, 0.16, 2, 10, 0, 1, 0.16, 2, -0.8)

% Computation of the Heston price using Milstein discretization and
% Monte-Carlo method

fprintf('Computation of the Heston price with the Milstein discretization + Monte-Carlo');
HestonPrice_MC(0, 1, 0.16, 2, 10, 0, 1, 0.16, 2, -0.8, 500, 1000)


With the close form formula I obtain : 0.0495 With the Monte-Carlo Simulation I obtain : 0.2300

If anyone of you have any idea about where this problem could come from, I would be very happy of it since I've spend 2 weeks trying to find the issue.

Thank you by advance ! :)

• Quick guess: your choice of parameters are "extreme". You have a vol of vol of 200%, since eta=2. The MC simulation can generate negative volatility which you compensate for. But if the negative volatility occurs frequently, then the error can grow quite large. – Olaf Jul 11 '15 at 17:42
• Another indicator is the Feller condition, given by $F = 2\lambda v_{\text{mean}}/\xi^2$. If $F>1$ then the volatility is strictly positive. But in your case $F=.08$. This is so small, that I suspect you are generating a lot of paths with negative volatility. – Olaf Jul 11 '15 at 17:45
• Thank you very much, I modifier my parameters and prices look closer. Just wondering : if I have a lot of negative volatility, I choose to set them to 0, so the volatility of the asset is lower, so the price show be lower in the case of the Monte-Carlo simulation, isn't it ? – Matt59 Jul 11 '15 at 21:42
• There are by now a huge number of papers on Heston and MC. The method you have chosen is not very good and there are many better methods. In particular, QE and my long-stepping method. See ssrn.com/abstract=1617187 . The code is downloadable from markjoshi.com – Mark Joshi Jul 12 '15 at 4:54
• @Matt59: actually, it's not that simple. The volatility is zero, but the stock price still grows at the mean rate. On top of that, volatility is negatively correlated with the stock price, so the paths for which negative volatility occurs are more likely to correspond to paths with relatively large stock prices. – Olaf Jul 12 '15 at 8:45