Calibrating and Pricing a Heston Stochastic-Local Volatility Model

# Mixed local-stochastic volatility model in Quantlib


# Mixed local-stochastic volatility model in Quantlib












At a conference the speaker mentioned that it is a standard approach today to use a mix of local and stochastic volatility model in equity, FX and interest rates. Can you please suggest the most intuitive and clear explanation of the process that goes from calibrating to market data to actual pricing of some exotic deal? Also, is there any library supporting this, like for example Quantlib? I cannot find any reference online.

## Answer by StackG (score 15)

https://quant.stackexchange.com/a/57147

Stochastic-Local Vol (SLV) is an attempt to mix the strengths and weaknesses of both Stochastic Vol and Local Vol models. Below, I'll quickly summarise each model and their strengths and weaknesses, and then discuss how SLV tries to improve things. Although there are many stochastic vol models, I limit the discussion here to the Heston model to keep things as short as possible. At the bottom, I've included some QuantLib-Python code that will calibrate, price options, and generate paths for exotic option pricing.

Local Vol

Local Vol typically refers to a generalisation of Black Scholes, where we assume a similar form of the underlying dynamics expect that a deterministic instantaneous volatility function is allowed to vary with both spot level $S$ and time $t$, so that risk-neutral dynamics obey \begin{align} dS = rS(t)dt + \sigma(S,t)S(t) dW_t \end{align}

This can correctly produce the prices of all observable vanilla options, if a continuous vol surface is observable (or can be interpolated) by setting \begin{align} \sigma(S,t) = \sqrt{{\frac {\frac {\partial C} {\partial T}} {{\frac 1 2} K^2 {} {\frac {\partial^2 C} {\partial K^2}}}}} \end{align}

The obvious strength of local vol is that it can exactly fit any observed vanilla surface, so once you've constructed the $\sigma(S,t)$ surface you don't need to worry about calibration or getting arbed.

However, the weakness is that it assumes deterministic volatility, so gets volatility dynamics badly wrong. This might not be a big problem for nearly-vanilla products like Asian options, but for products depending on vol-of-vol (eg. barrier options) or where forward vol is important (eg. forward-starting options), Local Vol produces prices that are far below what market participants charge.

Heston Stochastic Vol

The Heston model adds an additional stochastic driver for the instantaneous variance, so spot dynamics obey \begin{align} dS &= rS(t)dt + \sqrt{\nu(t)}S(t)dW^S_t \\ d\nu &= \kappa (\theta - \nu(t))dt + \epsilon \sqrt{\nu(t)}dW^{\nu}_t \end{align} and the two stochastic processes $dW^S_t, dW^\nu_t$ have correlation $\rho$

The variance equation is mean-reverting, so variance should move around the mean value of $\Theta$ with a reversion speed determined by $\kappa$. The 'vol-of-vol' term $\epsilon$ controls the amount of smile in the vol surface produced by this model, as it leads to increased vol when spot is already far from initial spot, and the correlation $\rho$ controls the skew of the surface.

This model does a better job of pricing vol-dependent options, but presents a problem of calibration. We only have 5 parameters that we can adjust, so we will not be able to produce a model that hits all of the vanilla option prices that are available. Instead, we calibrate the five parameters to observable prices using a least-squares procedure, chosing parameters that minimise the squared errors between the model proces and the observed prices (fortunately there is a semi-analytic pricing formula for vanilla options in Heston, so this is reasonably fast).

Stochastic-Local Vol

The insight in SLV is that we want to keep the dynamics from our stochastic vol model, but we need to adjust the average amount of volatility that the model picks up at each point on the $S, t$ surface so that it matches the amount of vol from the local vol model. This is achieved by adding a Leverage Function, $L(S,t)$ which scales up the vol that the stochastic vol model produces when it under-prices vanilla options, and scales down the the vol when it is too high (this is quite close to the ratio between the vol surfaces coming from the local vol model and the stochattic vol model, which is how I like to visualise it). In addition, a mixing fraction $\eta$ is usually added to calibrate between local and stochastic vol to price vol-dependent market exotics.

The resulting risk-neutral dynamics are \begin{align} dS &= rS(t)dt + \sqrt{\nu(t)}L(S,t)S(t)dW^S_t \\ d\nu &= \kappa (\theta - \nu(t))dt + \eta \epsilon \sqrt{\nu(t)}dW^{\nu}_t \end{align}

The procedure to calibrate this is:

- take an observed local vol surface and calculate the Dupire instantaneous vol

- calibrate a heston process to match it as well as you can

- then to pass both into the $L(S,t)$ calibration process - this is very complicated but fortunately we can delegate it to QuantLib

For sensible Heston parameters, this will give us back a model that exactly re-prices vanilla options, and the $\eta$ parameter can be adjusted from 0 to 1 to correctly price first-generation exotics (typically DNTs in FX, for example - https://arxiv.org/pdf/1911.00877.pdf)

QuantLib-Python Code

NOTE: all of this code requires the boilerplate code at the bottom of the post in order to run... but it's long, so I pushed it down

First, create and plot a vol surface using some random params in a heston process - let's assume this is the data the market has shown us today (imagine we don't know the process or the params that generated it...):

```
dates, strikes, vols, feller = create_vol_surface_mesh_from_heston_params(today, calendar, spot, 0.0225, 1.0, 0.0625, -0.25, 0.3, flat_ts, dividend_ts)

local_vol_surface = ql.BlackVarianceSurface(today, calendar, dates, strikes, vols, day_count)

# Plot the vol surface ...
plot_vol_surface(local_vol_surface)
```

Now, following the steps above:

- Calculate the Dupire instantaneous vol

```
spot_quote = ql.QuoteHandle(ql.SimpleQuote(spot))

local_vol_surface.setInterpolation("bicubic")
local_vol_handle = ql.BlackVolTermStructureHandle(local_vol_surface)
local_vol = ql.LocalVolSurface(local_vol_handle, flat_ts, dividend_ts, spot_quote)
local_vol.enableExtrapolation()

# Plot the Dupire surface ...
plot_vol_surface(local_vol, funct='localVol')
```

- Calibrate a Heston process (here, to make things interesting, let's assume we got the parameters slightly wrong... so the vol surface doesn't quite match)

```
# Create new heston model
v0 = 0.015; kappa = 2.0; theta = 0.065; rho = -0.3; sigma = 0.45; spot = 1007
feller = 2 * kappa * theta - sigma ** 2

heston_process = ql.HestonProcess(flat_ts, dividend_ts, spot_quote, v0, kappa, theta, sigma, rho)
heston_model = ql.HestonModel(heston_process)

# How does the vol surface look at the moment?
heston_handle = ql.HestonModelHandle(heston_model)
heston_vol_surface = ql.HestonBlackVolSurface(heston_handle)

# Plot the vol surface ...
plot_vol_surface([local_vol_surface, heston_vol_surface])
```

The heston surface is shown in orange - you can see it slightly misses the BS local vol surface, and the Leverage function will attempt to 'correct' for the difference

- Run the local vol fitting and calculate the leverage function

Calibrating the leverage function depends on a numerical accuracy parameter called `calibrationPaths`. The default is `2**15`, which runs quickly but produces quite a spikey leverage function. Increasing this parameter makes the leverage function converge to a smoother result, at the cost of increasing the runtime required (and I run out of memory above `2**19`), as shown here:

```
# Calibrate via Monte-Carlo
import time
end_date = ql.Date(1, 7, 2021)
generator_factory = ql.MTBrownianGeneratorFactory(43)

calibration_paths_vars = [2**15, 2**17, 2**19, 2**20]
time_steps_per_year, n_bins = 365, 201

for calibration_paths in calibration_paths_vars:
    print("Paths: {}".format(calibration_paths))
    stoch_local_mc_model = ql.HestonSLVMCModel(local_vol, heston_model, generator_factory, end_date, time_steps_per_year, n_bins, calibration_paths)

    a = time.time()
    leverage_functon = stoch_local_mc_model.leverageFunction()
    b = time.time()

    print("calibration took {0:2.1f} seconds".format(b-a))
    plot_vol_surface(leverage_functon, funct='localVol', plot_years=np.arange(0.1, 0.98, 0.1))
    plt.pause(0.05)
```

Now, let's create a path generator and generate paths from the Stoch Vol process:

```
num_paths = 25000
timestep = 32
length = 1
times = ql.TimeGrid(length, timestep)

stoch_local_process = ql.HestonSLVProcess(heston_process, leverage_functon)
dimension = stoch_local_process.factors()

rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(dimension * timestep, ql.UniformRandomGenerator()))
seq = ql.GaussianMultiPathGenerator(stoch_local_process, list(times), rng, False)

df_spot, df_vol = generate_multi_paths_df(seq, num_paths)

fig = plt.figure(figsize=(20,10))

plt.subplot(2, 2, 1)
plt.plot(df_spot.iloc[0:10].transpose())

plt.subplot(2, 2, 2)
plt.hist(df_spot[1.0])

plt.subplot(2, 2, 3)
plt.plot(df_vol.iloc[0:10].transpose())

plt.subplot(2, 2, 4)
plt.hist(df_vol[1.0])

plt.show()
```

- And at last, we can price options using these paths via Monte Carlo

Currently the QuantLib pricing engines for SLV haven't been reliably ported to Python, but I think they are coming soon!

```
# One year call at strike 100
(df_spot[1.0] - 100).clip_lower(0).mean()
```

- Actually one last thing...

You CAN actually price vanillas in QuantLib-Python using the analytical finite difference engine... and luckily you'll see this matches the MC price above closely:

```
slv_engine = ql.FdHestonVanillaEngine(heston_model, 400, 400, 100, 0, ql.FdmSchemeDesc.Hundsdorfer(), leverage_functon)

option.setPricingEngine(slv_engine)

option.NPV()
```

Boilerplate Code

```
import warnings
warnings.filterwarnings('ignore')

import QuantLib as ql
import numpy as np
import pandas as pd
import itertools

from scipy.stats import norm
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Some utility functions used later to plot 3D vol surfaces, generate paths, and generate vol surface from Heston params
def plot_vol_surface(vol_surface, plot_years=np.arange(0.1, 3, 0.1), plot_strikes=np.arange(70, 130, 1), funct='blackVol'):
    if type(vol_surface) != list:
        surfaces = [vol_surface]
    else:
        surfaces = vol_surface

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(plot_strikes, plot_years)

    for surface in surfaces:
        method_to_call = getattr(surface, funct)

        Z = np.array([method_to_call(float(y), float(x)) 
                      for xr, yr in zip(X, Y) 
                          for x, y in zip(xr,yr) ]
                     ).reshape(len(X), len(X[0]))

        surf = ax.plot_surface(X,Y,Z, rstride=1, cstride=1, linewidth=0.1)

    fig.colorbar(surf, shrink=0.5, aspect=5)

def generate_multi_paths_df(sequence, num_paths):
    spot_paths = []
    vol_paths = []

    for i in range(num_paths):
        sample_path = seq.next()
        values = sample_path.value()

        spot, vol = values

        spot_paths.append([x for x in spot])
        vol_paths.append([x for x in vol])

    df_spot = pd.DataFrame(spot_paths, columns=[spot.time(x) for x in range(len(spot))])
    df_vol = pd.DataFrame(vol_paths, columns=[spot.time(x) for x in range(len(spot))])

    return df_spot, df_vol

def create_vol_surface_mesh_from_heston_params(today, calendar, spot, v0, kappa, theta, rho, sigma, 
                                               rates_curve_handle, dividend_curve_handle,
                                               strikes = np.linspace(40, 200, 161), tenors = np.linspace(0.1, 3, 60)):
    quote = ql.QuoteHandle(ql.SimpleQuote(spot))

    heston_process = ql.HestonProcess(rates_curve_handle, dividend_curve_handle, quote, v0, kappa, theta, sigma, rho)
    heston_model = ql.HestonModel(heston_process)
    heston_handle = ql.HestonModelHandle(heston_model)
    heston_vol_surface = ql.HestonBlackVolSurface(heston_handle)

    data = []
    for strike in strikes:
        data.append([heston_vol_surface.blackVol(tenor, strike) for tenor in tenors])

    expiration_dates = [calendar.advance(today, ql.Period(int(365*t), ql.Days)) for t in tenors]
    implied_vols = ql.Matrix(data)
    feller = 2 * kappa * theta - sigma ** 2

    return expiration_dates, strikes, implied_vols, feller

# World State for Vanilla Pricing
spot = 100
rate = 0.0
today = ql.Date(1, 7, 2020)
calendar = ql.NullCalendar()
day_count = ql.Actual365Fixed()

# Set up the flat risk-free curves
riskFreeCurve = ql.FlatForward(today, rate, ql.Actual365Fixed())

flat_ts = ql.YieldTermStructureHandle(riskFreeCurve)
dividend_ts = ql.YieldTermStructureHandle(riskFreeCurve)
```