Local Stochastic Volatility with Monte-Carlo

I always imagined local stochastic volatility to be complicated, and thought it would be very slow to calibrate.

After reading a bit about it, I noticed that the calibration phase could just consist in calibrating independently a Dupire local volatility model and a stochastic volatility model the usual way.

One can then choose to compute on the fly the local volatility component (not equal the Dupire one, but including the stochastic adjustment) in the Monte-Carlo simulation to price a product.

There are two relatively simple algorithms that are remarkably similar, one by Guyon and Henry-Labordère in "The Smile Calibration Problem Solved":

And one from Van der Stoep, Grzelak & Oosterlee "The Heston Stochastic-Local Volatility Model: Efficient Monte Carlo Simulation":

In the particle method, the delta function is a regularizing kernel approximating the Dirac function. If we use the notation of the second paper, we have a = psi.

The methods are extremely similar, the evaluation of the expectation is slightly different, but even that is not very different. The disadvantage is that all paths are needed at each time step. As a payoff is evaluated over one full path, this is quite memory intensive.

Local Stochastic Volatility with Monte-Carlo

I always imagined local stochastic volatility to be complicated, and thought it would be very slow to calibrate.

After reading a bit about it, I noticed that the calibration phase could just consist in calibrating independently a Dupire local volatility model and a stochastic volatility model the usual way.

One can then choose to compute on the fly the local volatility component (not equal the Dupire one, but including the stochastic adjustment) in the Monte-Carlo simulation to price a product.

There are two relatively simple algorithms that are remarkably similar, one by Guyon and Henry-Labordère in "The Smile Calibration Problem Solved":

And one from Van der Stoep, Grzelak & Oosterlee "The Heston Stochastic-Local Volatility Model: Efficient Monte Carlo Simulation":

In the particle method, the delta function is a regularizing kernel approximating the Dirac function. If we use the notation of the second paper, we have a = psi.

The methods are extremely similar, the evaluation of the expectation is slightly different, but even that is not very different. The disadvantage is that all paths are needed at each time step. As a payoff is evaluated over one full path, this is quite memory intensive.

Heston, Schobel-Zhu, Bates, Double-Heston Fit

I did some experiments fitting Heston, Schobel-Zhu, Bates and Double-Heston to a real world equity index implied volatility surface. I used a global optimizer (differential evolution).

To my surprise, the Heston fit is quite good: the implied volatility error is less than 0.42% on average. Schobel-Zhu fit is also good (0.47% RMSE), but a bit worse than Heston. Bates improves quite a bit on Heston although it has 3 more parameters, we can see the fit is better for short maturities (0.33% RMSE). Double-Heston has the best fit overall but it is not that much better than simple Heston while it has twice the number of parameters, that is 10 (0.24 RMSE). Going beyond, for example Triple-Heston, does not improve anything, and the optimization becomes very challenging. With Double-Heston, I already noticed that kappa is very low (and theta high) for one of the processes, and kappa is very high (and theta low) for the other process: so much that I had to add a penalty to enforce constraints in my local optimizer. The best fit is at the boundary for kappa and theta. So double Heston already feels over-parameterized.

Heston volatility error

Schobel-Zhu volatility error

Bates volatility error

Double Heston volatility error

Another advantage of Heston is that one can find tricks to find a good initial guess for a local optimizer.

Update October 7: My initial fit relied only on differential evolution and was not the most stable as a result. Adding Levenberg-Marquardt at the end stabilizes the fit, and improves the fit a lot, especially for Bates and Double Heston. I updated the graphs and conclusions accordingly. Bates fit is not so bad at all.

Heston, Schobel-Zhu, Bates, Double-Heston Fit

I did some experiments fitting Heston, Schobel-Zhu, Bates and Double-Heston to a real world equity index implied volatility surface. I used a global optimizer (differential evolution).

To my surprise, the Heston fit is quite good: the implied volatility error is less than 0.42% on average. Schobel-Zhu fit is also good (0.47% RMSE), but a bit worse than Heston. Bates improves quite a bit on Heston although it has 3 more parameters, we can see the fit is better for short maturities (0.33% RMSE). Double-Heston has the best fit overall but it is not that much better than simple Heston while it has twice the number of parameters, that is 10 (0.24 RMSE). Going beyond, for example Triple-Heston, does not improve anything, and the optimization becomes very challenging. With Double-Heston, I already noticed that kappa is very low (and theta high) for one of the processes, and kappa is very high (and theta low) for the other process: so much that I had to add a penalty to enforce constraints in my local optimizer. The best fit is at the boundary for kappa and theta. So double Heston already feels over-parameterized.

Heston volatility error

Schobel-Zhu volatility error

Bates volatility error

Double Heston volatility error

Another advantage of Heston is that one can find tricks to find a good initial guess for a local optimizer.

Update October 7: My initial fit relied only on differential evolution and was not the most stable as a result. Adding Levenberg-Marquardt at the end stabilizes the fit, and improves the fit a lot, especially for Bates and Double Heston. I updated the graphs and conclusions accordingly. Bates fit is not so bad at all.

Second Cumulant of Heston

I recently stumbled upon an error in the various papers related to the Heston Cos method regarding the second cumulant. It is used to define the boundaries of the Cos method.

Letting phi be Heston characteristic function, the cumulant generating function is:
$$g(u) = \log(\phi(-iu))$$

And the second cumulant is defined a:
$$c_2 = g''(0)$$

Compared to a numerical implementation, the c_2 from the paper is really off in many use cases.

This is where Maxima comes useful, even if I had to simplify the results by hand. It leads to the following analytical formula:
$$c_2 = \frac{v_0}{4\kappa^3}\{ 4 \kappa^2 \left(1+(\rho\sigma t -1)e^{-\kappa t}\right)   + \kappa \left(4\rho\sigma(e^{-\kappa t}-1)-2\sigma^2 t e^{-\kappa t}\right)+\sigma^2(1-e^{-2\kappa t}) \}\\
+  \frac{\theta}{8\kappa^3} \{ 8 \kappa^3 t - 8 \kappa^2 \left(1+ \rho\sigma t + (\rho\sigma t-1)e^{-\kappa t}\right) + 2\kappa \left( (1+2e^{-\kappa t})\sigma^2 t+8(1-e^{-\kappa t})\rho\sigma \right) \\
+ \sigma^2(e^{-2\kappa t} + 4e^{-\kappa t}-5) \}$$

In contrast, the paper formula was:

I saw this while trying to calibrate Heston on a bumped surface: the results were very different with the Cos method than with the other methods. The short maturities were mispriced, except if one pushed the truncation level L to 24 (instead of the usual 12), and as a result one would also need to significantly raise the number of points used in the Cos method. With the corrected formula, it works well with L=12.

Here is an example of failure on a call option of strike, spot 1.0 and maturity 1.0 and not-so-realistic Heston parameters $$\kappa=0.1, \theta=1.12, \sigma=1.0, v_0=0.2, \rho=-0.377836$$ using 200 points:

New formula: 0.1640581405
Paper formula: 0.1743425406
N=30 with 10000 points: 0.1640581423

Update March 2014 - this is now described in my paper Fourier Integration and Stochastic Volatility Calibration.

Second Cumulant of Heston

I recently stumbled upon an error in the various papers related to the Heston Cos method regarding the second cumulant. It is used to define the boundaries of the Cos method.

Letting phi be Heston characteristic function, the cumulant generating function is:
$$g(u) = \log(\phi(-iu))$$

And the second cumulant is defined a:
$$c_2 = g''(0)$$

Compared to a numerical implementation, the c_2 from the paper is really off in many use cases.

This is where Maxima comes useful, even if I had to simplify the results by hand. It leads to the following analytical formula:
$$c_2 = \frac{v_0}{4\kappa^3}\{ 4 \kappa^2 \left(1+(\rho\sigma t -1)e^{-\kappa t}\right)   + \kappa \left(4\rho\sigma(e^{-\kappa t}-1)-2\sigma^2 t e^{-\kappa t}\right)+\sigma^2(1-e^{-2\kappa t}) \}\\
+  \frac{\theta}{8\kappa^3} \{ 8 \kappa^3 t - 8 \kappa^2 \left(1+ \rho\sigma t + (\rho\sigma t-1)e^{-\kappa t}\right) + 2\kappa \left( (1+2e^{-\kappa t})\sigma^2 t+8(1-e^{-\kappa t})\rho\sigma \right) \\
+ \sigma^2(e^{-2\kappa t} + 4e^{-\kappa t}-5) \}$$

In contrast, the paper formula was:

I saw this while trying to calibrate Heston on a bumped surface: the results were very different with the Cos method than with the other methods. The short maturities were mispriced, except if one pushed the truncation level L to 24 (instead of the usual 12), and as a result one would also need to significantly raise the number of points used in the Cos method. With the corrected formula, it works well with L=12.

Here is an example of failure on a call option of strike, spot 1.0 and maturity 1.0 and not-so-realistic Heston parameters $$\kappa=0.1, \theta=1.12, \sigma=1.0, v_0=0.2, \rho=-0.377836$$ using 200 points:

New formula: 0.1640581405
Paper formula: 0.1743425406
N=30 with 10000 points: 0.1640581423

Update March 2014 - this is now described in my paper Fourier Integration and Stochastic Volatility Calibration.

Maxima for Symbolic Calculus

A few years ago, I found an interesting open source symbolic calculus software called Xcas. It can however be quickly limited, for example, it does not seem to work well to compute Taylor expansions with several embedded functions.
Google pointed me to another popular open source package, Maxima. It looks a bit rudimentary (command like interface), but formulas can actually be very easily exported to latex with the tex command. Here is a simple example:

(%i14) D(x):=sqrt((lambda-rho*eta*x)^2+(-x^2+x)*eta^2);
                                              2       2         2
(%o14)       D(x) := sqrt((lambda - rho eta x)  + (- x  + x) eta )

(%i15) G(x) := (lambda - rho*eta*x - D(x))/(lambda - rho*eta*x +D(x));
                               lambda - rho eta x - D(x)
(%o15)                 G(x) := -------------------------
                               lambda - rho eta x + D(x)

(%i16) tex(taylor((1-exp(-t*D(x)))/(1-G(x)*exp(-t*D(x)))*(lambda - rho*eta*x - D(x)),x,0,3));
$$-{{\left(e^{t\,\lambda}-1\right)\,\eta^2\,x}\over{2\,e^{t\,\lambda}
 \,\lambda}}+{{\left(\left(4\,e^{t\,\lambda}\,t\,\eta^3\,\rho+\left(4
 \,\left(e^{t\,\lambda}\right)^2-4\,e^{t\,\lambda}\right)\,\eta^2
 \right)\,\lambda^2+\left(\left(-4\,\left(e^{t\,\lambda}\right)^2+4\,
 e^{t\,\lambda}\right)\,\eta^3\,\rho-2\,e^{t\,\lambda}\,t\,\eta^4
 \right)\,\lambda+\left(\left(e^{t\,\lambda}\right)^2-1\right)\,\eta^
 4\right)\,x^2}\over{8\,\left(e^{t\,\lambda}\right)^2\,\lambda^3}}+{{
 \left(\left(8\,\left(e^{t\,\lambda}\right)^2\,t^2\,\eta^4\,\rho^2-16
 \,\left(e^{t\,\lambda}\right)^2\,t\,\eta^3\,\rho\right)\,\lambda^4+
 \left(16\,\left(e^{t\,\lambda}\right)^2\,t\,\eta^4\,\rho^2+\left(-8
 \,\left(e^{t\,\lambda}\right)^2\,t^2\,\eta^5+\left(16\,\left(e^{t\,
 \lambda}\right)^3-16\,\left(e^{t\,\lambda}\right)^2\right)\,\eta^3
 \right)\,\rho+16\,\left(e^{t\,\lambda}\right)^2\,t\,\eta^4\right)\,
 \lambda^3+\left(\left(-16\,\left(e^{t\,\lambda}\right)^3+16\,\left(e
 ^{t\,\lambda}\right)^2\right)\,\eta^4\,\rho^2+\left(-16\,\left(e^{t
 \,\lambda}\right)^2-8\,e^{t\,\lambda}\right)\,t\,\eta^5\,\rho+2\,
 \left(e^{t\,\lambda}\right)^2\,t^2\,\eta^6+\left(-8\,\left(e^{t\,
 \lambda}\right)^3+8\,e^{t\,\lambda}\right)\,\eta^4\right)\,\lambda^2
 +\left(\left(12\,\left(e^{t\,\lambda}\right)^3-12\,e^{t\,\lambda}
 \right)\,\eta^5\,\rho+\left(2\,\left(e^{t\,\lambda}\right)^2+4\,e^{t
 \,\lambda}\right)\,t\,\eta^6\right)\,\lambda+\left(-2\,\left(e^{t\,
 \lambda}\right)^3-\left(e^{t\,\lambda}\right)^2+2\,e^{t\,\lambda}+1
 \right)\,\eta^6\right)\,x^3}\over{32\,\left(e^{t\,\lambda}\right)^3
 \,\lambda^5}}+\cdots $$

Regarding Taylor expansion, there seems to be quite a few options possible, but I found that the default expansion was already relatively easy to read. XCas produced less readable expansions, or just failed.

Maxima for Symbolic Calculus

A few years ago, I found an interesting open source symbolic calculus software called Xcas. It can however be quickly limited, for example, it does not seem to work well to compute Taylor expansions with several embedded functions.
Google pointed me to another popular open source package, Maxima. It looks a bit rudimentary (command like interface), but formulas can actually be very easily exported to latex with the tex command. Here is a simple example:

(%i14) D(x):=sqrt((lambda-rho*eta*x)^2+(-x^2+x)*eta^2);
                                              2       2         2
(%o14)       D(x) := sqrt((lambda - rho eta x)  + (- x  + x) eta )

(%i15) G(x) := (lambda - rho*eta*x - D(x))/(lambda - rho*eta*x +D(x));
                               lambda - rho eta x - D(x)
(%o15)                 G(x) := -------------------------
                               lambda - rho eta x + D(x)

(%i16) tex(taylor((1-exp(-t*D(x)))/(1-G(x)*exp(-t*D(x)))*(lambda - rho*eta*x - D(x)),x,0,3));
$$-{{\left(e^{t\,\lambda}-1\right)\,\eta^2\,x}\over{2\,e^{t\,\lambda}
 \,\lambda}}+{{\left(\left(4\,e^{t\,\lambda}\,t\,\eta^3\,\rho+\left(4
 \,\left(e^{t\,\lambda}\right)^2-4\,e^{t\,\lambda}\right)\,\eta^2
 \right)\,\lambda^2+\left(\left(-4\,\left(e^{t\,\lambda}\right)^2+4\,
 e^{t\,\lambda}\right)\,\eta^3\,\rho-2\,e^{t\,\lambda}\,t\,\eta^4
 \right)\,\lambda+\left(\left(e^{t\,\lambda}\right)^2-1\right)\,\eta^
 4\right)\,x^2}\over{8\,\left(e^{t\,\lambda}\right)^2\,\lambda^3}}+{{
 \left(\left(8\,\left(e^{t\,\lambda}\right)^2\,t^2\,\eta^4\,\rho^2-16
 \,\left(e^{t\,\lambda}\right)^2\,t\,\eta^3\,\rho\right)\,\lambda^4+
 \left(16\,\left(e^{t\,\lambda}\right)^2\,t\,\eta^4\,\rho^2+\left(-8
 \,\left(e^{t\,\lambda}\right)^2\,t^2\,\eta^5+\left(16\,\left(e^{t\,
 \lambda}\right)^3-16\,\left(e^{t\,\lambda}\right)^2\right)\,\eta^3
 \right)\,\rho+16\,\left(e^{t\,\lambda}\right)^2\,t\,\eta^4\right)\,
 \lambda^3+\left(\left(-16\,\left(e^{t\,\lambda}\right)^3+16\,\left(e
 ^{t\,\lambda}\right)^2\right)\,\eta^4\,\rho^2+\left(-16\,\left(e^{t
 \,\lambda}\right)^2-8\,e^{t\,\lambda}\right)\,t\,\eta^5\,\rho+2\,
 \left(e^{t\,\lambda}\right)^2\,t^2\,\eta^6+\left(-8\,\left(e^{t\,
 \lambda}\right)^3+8\,e^{t\,\lambda}\right)\,\eta^4\right)\,\lambda^2
 +\left(\left(12\,\left(e^{t\,\lambda}\right)^3-12\,e^{t\,\lambda}
 \right)\,\eta^5\,\rho+\left(2\,\left(e^{t\,\lambda}\right)^2+4\,e^{t
 \,\lambda}\right)\,t\,\eta^6\right)\,\lambda+\left(-2\,\left(e^{t\,
 \lambda}\right)^3-\left(e^{t\,\lambda}\right)^2+2\,e^{t\,\lambda}+1
 \right)\,\eta^6\right)\,x^3}\over{32\,\left(e^{t\,\lambda}\right)^3
 \,\lambda^5}}+\cdots $$

Regarding Taylor expansion, there seems to be quite a few options possible, but I found that the default expansion was already relatively easy to read. XCas produced less readable expansions, or just failed.