Comparing Heston Closed-Form Pricing with Milstein Monte Carlo

# Simulation of Heston process


# 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 ! :)