# Euler discretization of Heston SDE in Mathematica

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]]


change the discretization and use the QE-M approach: Andersen (2006) the bias is way smaller than the one of the simple Euler. further u can try to use control variates/anthitetic numbers to reduce the sample variance.

Some simple improvements:

1) Replace the Euler discretization approximation of the volatility to a Milstein discretization approximation. See e.g. these notes by Rouah.

2) 100 Paths is a very low number of paths, and leads to a big standard error in your estimate. So this should be increased by a factor of ~100.

3) You should use some form of variance reduction. Antithetic variables are easy to implement and give a great improvement on your standard error.

As Phun mentioned, there are a number of more complicated approaches that you could use to construct your paths. These reduce your bias, since these can approaches can typically avoid the variance from becoming negative.

But I would suggest trying the above simplifications first, because the bias should not be extremely big. With this scheme you should be able to get pretty close to the implied volatilities as determined by the Fourier pricing method.

And finally, some sanity checks: did you calibrate the Heston model using a zero rate as well? Do you exponentiate the final HestonPaths to get the "actual" Heston Paths? (since you are approximating $d\log[S]$).

Thanks for the responses. I'm still puzzling over this.

Here is an implementation which uses both Milstein discretization as well as antithetic variables.

The code constructs the volatility skew for a T = 574 day call with initial forward price G0 = 570.856 and rate of interest r = 0.05327. The volatility skew is clearly crazy: it is concave with a positive skew whereas the Heston model should give a convex smile and a negative skew.

I've set the number of increments to be n = 100 and NPaths = 100 just to demonstrate the code. Increasing to n = 1000 increments with NPaths = 10000 takes about an hour to execute and does not seem to improve the result.

\[Rho] = -0.4042;
v0 = 0.2577^2;
\[Kappa] = 0.2656;
vbar = .1851;
\[Sigma]v = 0.2992;
n = 100;
NPaths = 100;

HestonPaths[G0_, T_] :=
Module[{XPaths, X1, v1, X2, v2, dt, dw},
dt = T/n;
dw = RandomVariate[
BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]], {NPaths, n}];
XPaths = {};
Do[
X1 = {Log[G0]};
v1 = {v0};
X2 = {Log[G0]};
v2 = {v0};
Do[
v1 = Append[v1,
Abs[\[Kappa] (vbar - Last[v1]) dt -
1/4 \[Sigma]v^2 dt + (Sqrt[Last[v1]] +
1/2 \[Sigma]v dw[[idx, i]][[1]])^2]];
X1 = Append[X1,
Last[X1] - 1/2 Last[v1] dt + Sqrt[Last[v1]] dw[[idx, i]][[2]]];
v2 = Append[v2,
Abs[\[Kappa] (vbar - Last[v2]) dt -
1/4 \[Sigma]v^2 dt + (Sqrt[Last[v1]] -
1/2 \[Sigma]v dw[[idx, i]][[1]])^2]];
X2 = Append[X2,
Last[X2] - 1/2 Last[v2] dt - Sqrt[Last[v2]] dw[[idx, i]][[2]]];
, {i, 1, n}];
XPaths = Append[XPaths, X1];
XPaths = Append[XPaths, X2];
, {idx, 1, NPaths}];
Exp[XPaths]
]

T = 574/365;
G0 = 570.856;
r = 0.05327;

Paths = HestonPaths[G0, T];

CallPrice[K_] :=
Block[{CF, \[CurlyPhi]VanillaCall, G},
G = Table[Part[Paths, i, n], {i, 1, NPaths}];
\[CurlyPhi]VanillaCall[S_] := Max[0, S - K];
CF = Map[\[CurlyPhi]VanillaCall, G];
Exp[-r T] Mean[CF]
]

BCallIV[G_, K_, r_, T_, value_] :=
FinancialDerivative[{"European", "Call"}, {"StrikePrice" -> K,
"Expiration" -> T, "Value" -> value}, {"InterestRate" -> r,
"CurrentPrice" -> G, "Dividend" -> r}, "ImpliedVolatility"
];

ListPlot[Table[{k, BCallIV[G0, k, r, T, CallPrice[k]]}, {k, 300, 700,
10}]]

• Tricky to debug, because I dont have Mathematica. Is writing T as a fractional a potential cause for slowing down your code? I know mathematica likes to keep fractional expressions for as long as possible. 1 hour is quite excessive. Also, shouldn't the Dividend be set to zero in BCallIV? – Olaf Apr 3 '15 at 10:04
• Another thing you could check is the distribution of the end-points of the generated Paths (e.g. a histogram). This should resemble something log-normal like, but with fatter tails. – Olaf Apr 3 '15 at 10:09
• I'm modeling forward prices $G(t,T)$ rather than stock prices $S(t)$ which is why I have set $\delta = r$. The problem most definitely lies in the function HestonPaths because I have tested the same code using an implementation of the CEV model and everything works as expected. – user11881 Apr 3 '15 at 19:10

I seem to get much better results if I replace

BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho]]


with

BinormalDistribution[{Sqrt[dt], Sqrt[dt]}, \[Rho] dt]


This is very strange because the Mathematica documentation clearly states that BinormalDistribution takes the correlation $\rho$ as input, whereas the covariance is given by $\rho \, dt$.

I will try running this modified code with $10^4$ sample paths and $10^4$ time increments.

Is this overkill for getting decent statistics?