Euler discretization of Heston SDE in Mathematica

Question

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

Answer 1

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.

Answer 2

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]$).

Answer 3

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

Answer 4

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?