Below is an implementation of the numerical solution of the Heston SDE using Euler discretization. It takes under a second to run on Mathematica.
The calibration parameters give a good fit to the volatility surface using the characteristic function/Fourier transform technique.
I am trying to use the below code to price an exotic derivative by MC simulation but I am unable to match the volatility surface as a first step. I suspect the code is simply taking too long to converge.
Are there any quick fixes that can speed this code up significantly?
\[Rho] = -0.4042;
v0 = 0.2577^2;
\[Kappa] = 0.2656;
vbar = .1851;
\[Sigma]v = 0.2992;
n = 100;
NPaths = 100;
Tmax = 574/365;
dt = Tmax/n;
dw = RandomVariate[
BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]], {NPaths, n}];
HestonPaths[G0_] :=
Module[{XPaths, X, v},
XPaths = {};
Do[
X = {Log[G0]};
v = {v0};
Do[
v = Append[v,
Abs[Last[v] + \[Kappa] (vbar - Last[v]) dt +
Sqrt[Last[v]] \[Sigma]v dw[[idx, i]][[1]]]];
X = Append[X,
Last[X] - 1/2 Last[v] dt + Sqrt[Last[v]] dw[[idx, i]][[2]]];
, {i, 1, n}];
XPaths = Append[XPaths, X];
, {idx, 1, NPaths}];
Exp[XPaths]
]
ListLinePlot[HestonPaths[500]]