Monte-Carlo Calibration of the Heston Stochastic Local Volatiltiy Model

January 10, 2016September 14, 2017 hpcquantlib13 Comments

Solving the Fokker-Planck equation via finite difference methods is not the only way to calibrate the Heston stochastic local volatility model

$\begin{array}{rcl} d \ln{S_t} &=& \left(r_t - q_t - \frac{1}{2}L\left(S_t, t\right)^2\nu_t \right)dt + L(S_t, t)\sqrt{\nu_t} dW_t^S \\ \nonumber d\nu_t &=& \kappa\left(\theta - \nu_t\right)dt + \eta\sigma\sqrt{\nu_t}dW_t^\nu \\ \nonumber \rho dt &=& dW_t^\nu dW_t^S \end{array}$

The basic equation to calibrate the leverage function $L(S_t,t)$ for a local volatility surface $\sigma_{LV}(S_t,t)$ and a set of Heston parameters $\{\nu_0, \kappa, \theta, \sigma, \rho, \eta \}$ is given by

$L(S_t, t)=\frac{\sigma_{LV}(S_t, t)}{\sqrt{\mathop{\mathbb{E}}[\nu_t|S=S_t]}}= \sigma_{LV}(S_t, t)\sqrt{\frac{\int_{\mathop{\mathbb{R}}^+} p(S_t,\nu, t) d\nu}{\int_{\mathop{\mathbb{R}}^+}\nu p(S_t,\nu,t)d\nu}}$

Key problem here is to calculate the expectation value $\mathop{\mathbb{E}}[\nu_t|S=S_t]$ . This can either be done via the Fokker-Planck equation as outlined in [3] and the references in there or via Monte-Carlo simulations as shown in [2]. Getting the solution of the Fokker-Planck equation via finite difference methods is anything but easy if the Feller constraint is violated. In this case the Monte-Carlo approach promises less numerical problems to overcome. Starting point for an efficient Monte-Carlo calibration is a fast and accurate simulation scheme for a stochastic local volatility (SLV) model. The variance part of the SLV can be sampled exactly using the non-central $\chi^2$ distribution.

$Pr(\nu_{t+\Delta t} < \nu | \nu_t) = F_{\chi^{\prime }2}\left(\nu \frac{4\kappa}{\eta^2\sigma^2\left(1-e^{-\kappa \Delta t}\right)}; \frac{4\kappa\theta}{\eta^2\sigma^2}, \frac{4\kappa}{\eta^2\sigma^2\left(1-e^{-\kappa \Delta t}\right)}\nu_t e^{-\kappa \Delta t}\right)$

The boost library supports the inverse of the cumulated non-central $\chi^2$ distribution but the calculation itself is rather slow. The Heston QE scheme outlined in [1] comes with an fast and accurate sampling scheme for the square root process. In the case $\gamma_1=\gamma_2=\frac{1}{2}$ it reads

$\begin{array}{rcl} m &=& \theta+(\nu_t - \theta)e^{-\kappa \Delta t} \\[6pt] \nonumber \psi &=& \frac{\nu_t \frac{\eta^2\sigma^2 e^{-\kappa\Delta t}}{\kappa}\left(1-e^{-\kappa\Delta t} \right) + \frac{\theta\eta^2\sigma^2}{2\kappa}\left(1-e^{-\kappa\Delta t}\right)^2}{m^2}\\[6pt] \nonumber \beta &=& \frac{2}{\psi} - 1 + \sqrt{\frac{2}{\psi}\left(\frac{2}{\psi}-1 \right)} \\[6pt] \nonumber p &=& \frac{\psi-1}{\psi+1} \\[6pt] \nonumber Z_\nu &\sim& \mathcal{N}(0,1) , u \sim \mathcal{U}(0,1) \\[6pt] \nonumber \nu_{t+\Delta t} &=& \begin{cases} \frac{m}{1+\beta}\left(\sqrt{\beta}+Z_\nu\right)^2 & \text{if } \psi < 1.5 \\[6pt] \begin{cases} 0, & \text{if } u \leq p \\[6pt] \ln{\left(\frac{1-p}{1-u}\right)}\frac{m}{1-p} & \text{otherwise} \end{cases} & \text{if } \psi \geq 1.5 \end{cases} \end{array}$

The authors in [2] suggesting a sampling scheme for $X_t = \ln{S_t}$ . The scheme below is a bit closer to the original QE scheme and without a rigorous proof it has shown slightly smaller discretisation errors for the examples.

$\begin{array}{rcl} X_{t+\Delta t} &=& X_t + \left(r_t - q_t - \frac{1}{4}\left(\nu_t+\nu_{t+\Delta t}\right)L(S_t, t)^2 \right)\Delta t \\[6pt] \nonumber &&+ \frac{\rho}{\sigma}L(S_t, t) \left(\nu_{t+\Delta t} -\nu_t - \kappa\theta\Delta t + \frac{1}{2}\left(\nu_t + \nu_{t + \Delta t}\right)\kappa\Delta t\right) \\[6pt] \nonumber &&+\sqrt{\frac{1}{2}\left(1-\rho^2\right)\left(\nu_t + \nu_{t+\Delta t}\right) Z_x \Delta t} \\[6pt] \nonumber Z_x &\sim& \mathcal{N}(0,1) \end{array}$

To get from the Monte-Carlo paths to the leverage function the authors in [2] using a binning technique in every time step. QuantLib supports two ways to generate the random draws, either based on i.i.d normal variables or based on Sobol quasi Monte-Carlo simulations with brownian bridges.

Test case definition: Feller constraint is violated.

$\begin{array}{rcl} S_{0} &=& 100 \\ \nonumber r_t &=& 5\% \\ \nonumber q_t &=& 2\% \\ \nonumber \sigma_{LV}(S_t, t) &=& 30\% \\ \nonumber \kappa &=& 1.0\\ \nonumber \theta &=& 0.06\\ \nonumber \rho &=& -0.75\\ \nonumber \sigma &=& 0.4\\ \nonumber \nu_{0} &=& 0.09 \end{array}$

Effectively the leverage function $L(S_t, t)$ will remove the skew and term structure introduced by the Heston model and shifts volatility to 30%. In order to measure the calibration quality the price mismatches expressed in volatility basis point differences are calculated for OTM calls and puts with prices above 0.02 and for maturity from one month until two year. The SLV prices are calculated by solving the Feynman-Kac backward equation using finite difference methods.

avgpriceerror

The diagram above shows the average pricing error for different number of Monte-Carlo calibration paths and for different bin sizes. As expected the optimal bin size depends on the number of Monte-Carlo simulations. Especially deep OTM options benefit from a large number of bins. The usage of Sobol quasi Monte-Carlo simulations with brownian bridges does not give a significant advantage.

The average pricing error is decreasing with the square root of the number of Monte-Carlo calibration paths if one always choses the optimal bin size as shown in the diagram below. Quasi Mont-Carlo path generation with brownian bridges gives a small advanage.

optimalbinsize

As expected the leverage function itself exposes some Monte-Carlo noise, see diagram below. The corresponding leverage function based on the finite difference method is smooth but the overall calibration error is of the same size.

mc_leverage

The calibration for different mixing factors $\eta$ is straight forward. The diagram below e.g. shows the prices of a double no touch option for different distances of the (symemtric) lower and upper barriers expressed in terms of the corresponding Black-Scholes prices. This diagram type is derived from figure 8.8 in [4].

moustache

The source code is available on github as part of the test suite HestonSLVModelTest.

[1] Leif Andersen, Efficient Simulation of the Heston Stochastic Volatility Model

[2] Anthonie Van der Stoep, Lech Grzelak, Cornelis Oosterlee,
The Heston Stochastic-Local Volatility Model: Efficient Monte Carlo Simulation

[3] Johannes Goettker-Schnetmann, Klaus Spanderen, Calibrating the Heston Stochastic Local Volatility Model using the Fokker-Planck Equation

[4] Iain J. Clark, Foreign Exchange Option Pricing: A Practitioner’s Guide

13 thoughts on “Monte-Carlo Calibration of the Heston Stochastic Local Volatiltiy Model”

Philipp says:

January 14, 2016 at 9:19 pm

Great! Thank you very much. I am asking one more time regarding forward smiles LSV vs LV. Due to the Mimicking Theorem shouldn’t all future transition probabilities be equal in LV and LSV model? If all transition probabilities (including future ones) are equal for LV and LSV model, shouldn’t then forward smiles be equal of LV and LSV Independent of mixing fraction? I guess I am still missing the point.

Reply
Philipp says:

January 14, 2016 at 9:20 pm

Transition probabilities of LV and LSV are of course only equal when the varianse is integrated out for the LSV one.

Reply
hpcquantlib says:

January 15, 2016 at 5:06 pm

I think the future transition probabilities \gamma(x_1, t_1, x_2, t_2) are not necessarily the same even if the terminal distribution \gamma(x_1=x_0, t_1=0, x, t) is given by the implied volatility surface, which is reproduced by both models. Different propagators \gamma(x_1, t_1, x_2, t_2) lead to different forward skews and different prices.

Reply
Philipp says:

January 17, 2016 at 12:11 pm

The terminal distribution \gamma(x_1=x_0, t_1=0, x, t) is equal of LV and LSV model at all times, not only at the times where vanilla prices are given, right? For each pair, (S_t, t), the probability density of LV model: p(S_t,t) equals the one of LSV model :Integral p(S_t,t,v_t) dv_t. If this is true for all (S_t, t), I wonder how the future transition probabilities can be different. In my MC simulations, LV and LSV Forward smiles conditional on S_t = x are also equal. Is my mistake maybe that I interpolate the local vol Surface for each Point (S_t,t) to calculate the local volatility correction at each pont L(S_t,t) no matter which Point S(S_t,t) is given in the MC simulation?

Reply
hpcquantlib says:

January 17, 2016 at 9:13 pm

I fully agree, the terminal distribution \gamma(x_1=x_0, t_1=0, x, t) is equal for LV and LSV for all {x, t} . But from this IMO we can not infer that the propagator \gamma(x_1, t_1, x_2, t_2) is the same for all {x_1, t_1, x_2, t_2}.

Based on the test case described above I’ve written a small Monte-Carlo simulation to calculate the forward implied volatility for an option starting in one year from today and with maturity in two years time. The source code is available on github, repository https://github.com/jschnetm/quantlib/tree/slv, test case HestonSLVModel::testForwardSkew(). Even though the local volatility is constant everywhere (30%) the LSV model exhibits a skew for forward starting options if \eta is greater than zero. Please find the corresponding diagram here https://hpcquantlib.files.wordpress.com/2016/01/forwardskew.png .

Shall we compare the leverage functions e.g. for the test case in the blog entry in order to find the differences in the implementations?

Reply
Philipp says:

January 18, 2016 at 11:26 pm

You are right, I found a mistake and now get similar results like you. It also makes sense to me now from a conceptional point of view, I believe. I would appreciate if you let me know if the following makes sense to you: I analyse forward smiles, conditional on spot being at a certain level at the forward time T1 for options expiring at T2. So I don’t look at average smiles for certain level of moneyness like you did. For the SV model, the uncertainty about the instantaneous volatility at T1 increases with higher forward times T1. As ITM and OTM options have positive vomma, they are convext in implied volatility which is an approximation of avarage instantaneous volatility. Jensen’s inequality -> ITM and OTM Options have higher implied volatility than the average of instantaneous volatility. This makes the forward smile for SV models more convex the higher the forward time T1, as the uncertainty about instantaneous volatility increases. For the SLV model, the expectation of instantaneous variance is equal to the LV model. However, there is still stochasticity in the instantaneous variance. This makes Jensen’s inequality work again, and thus adds convexity to the Forward smiles the higher the forward time T1. What I don’t understand then is however, why you get lower implied vols for higher strikes and higher implied vols for lower strike in your graph https://hpcquantlib.files.wordpress.com/2016/01/forwardskew.png . Shouldn’t the smile be rather symmetric? Thank you very much for giving me insides during the last few weeks, you helped me very much in understanding these topics.

Reply
hpcquantlib says:

January 25, 2016 at 4:34 pm

I think the shape of the skew/smile of the LSV model is given by the Heston model parameters, e.g. if I use different correlations for the Heston model I get a forward smile (\rho=0.0) or left/right forward skew, please see https://hpcquantlib.files.wordpress.com/2016/01/forwardskewrho.png. Unfortunately I don’t have a short descriptive explanation for it other than that the local volatility forward skew flattens out and the forward skew/smile is dominated by the underlying Heston model parameters after some time.

Reply
Philipp says:

January 27, 2016 at 12:11 pm

Okay, thank you very much

Reply
Philipp says:

March 10, 2016 at 5:09 pm

I think I finally found an intuitive reasoning for it. The “Mimicking” or Gyöngy’s Theorem only holds for the expectation at t = 0 but not for the conditional expectation given a certain spot level in the future, i.e. not on each simulated path. This mean it holds for the expectation of variance at a certain point S_t and t over all possible paths but does not hold for any given path as soon as t > 0 and conditional on the spot level given by that path. For non path dependent options, this is sufficient to result in equal smiles, but not for path dependent options. This is because path dependent options have an expected payoff equal to the average of payoffs resulting from each path, and not the payoff resulting from the path that used the average variance.

Reply
vgdev says:

June 14, 2017 at 12:31 pm

Where can I find the code for this?

Reply
- hpcquantlib says:
  
  June 15, 2017 at 7:53 am
  
  Code is in QuantLib since release 1.8. Good starting points are the test cases HestonSLVModelTest::testMonteCarloCalibration and HestonSLVModelTest::testMoustacheGraph.
  
  Reply
Manoj Kumar Agrawal says:

September 11, 2017 at 2:05 pm

Has the source code been removed from Github for HestonSLVModelTest?

Reply
hpcquantlib says:

September 14, 2017 at 5:30 am

All code has been moved to the official QuantLib repository, please see
https://github.com/lballabio/QuantLib/blob/master/test-suite/hestonslvmodel.cpp

Thanks for the hint, I’ve change the link in the blog entry.

Reply