Computing Skew Stickiness Ratios in the Heston Model

# Calculating the skew-stickiness ratio in Affine Forward Variance Models


# Calculating the skew-stickiness ratio in Affine Forward Variance Models












I am trying to calculate the skew-stickiness ratio implied by a number of different stochastic volatility models. However, I am unsure about the right approach, as there seems to be almost no literature about this topic, and certainly no public code (except for the quadratic Heston, as can be found here). A recent paper details two equations which seem doable to implement myself in python. Focusing on Eq. 29, given by

\begin{equation} R_t^{N}(\tau) = \frac{1}{S_t(\tau)} \lim_{h\downarrow 0}\frac{1}{h} \left( \sigma^{\tau,N}\!\left( t,\ \Big\{\, \xi_{t}^{u_i} + h\,\frac{\lambda\!\big(t,u_i,\xi_t^{u_i}\big)\,\rho_{S\xi}(t)} {\sqrt{\xi_{t}^{t}}} \,\Big\}_{\,t<u_i<t+\tau} \right) - \sigma^{\tau,N}\!\left( t,\ \{\xi_{t}^{u_i}\}_{\,t<u_i<t+\tau} \right) \right). \end{equation} Which is also equal to Eq. 42 of the same paper (for the rough Heston) when using $H=0.5$, given by

\begin{equation} \mathcal{R}^{\mathrm{RH}}_t(\tau)\;\approx\; \frac{1}{S^{\mathrm{RH}}_t(\tau)}\,\frac{1}{h}\! \left( \sigma^{\tau,N}\!\left( t,\ \Big\{\, \xi^{u_i}_t + h\,\eta\,\rho\,(u_i-t)^{H-\tfrac{1}{2}}\, f_H\!\big(-\kappa\,(u_i-t)^{H+\tfrac{1}{2}}\big) \,\Big\}_{\,t<u_i<t+\tau} \right) - \sigma^{\tau,N}\!\left( t,\ \{\xi^{u_i}_t\}_{\,t<u_i<t+\tau} \right) \right). \end{equation} With $f_{\frac{1}{2}}(x) = e^x$

If my understanding of the theory is correct, we have for the Heston model that $\xi_0(u)$ is given by

\begin{equation} \xi_0(u) = v_0 e^{-\kappa u} + \theta\!\left(1 - e^{-\kappa u}\right) = \theta + (v_0 - \theta)\,e^{-\kappa u}. \end{equation} So, in the case of $u=0$ $\xi_u(u)$ reduces to $v_0$. then, if I correctly follow the paper, I get $\lambda(...) = h\times {\rho \times \sigma \times \exp(-\kappa u)} $ and $h = 0.01 \times \sqrt{T}$ with $T$ time to maturity of the option.

Then, we have the following characteristic function for the heston

```
def forward_heston_char(phi, v0, T,
                        kappa, theta, sigma, rho, xi0 = None,
                        n_quad=64):
    """
    Vectorized AFV-form Heston CF with
      g0(u)=xi0(u),
      f2(s)=κθ,
      ψ2' = -κψ2+F(iφ,ψ2), ψ2(0)=0.
    Accepts scalar or array phi.
    
    Returns array of E[e^{iφ ln S_T}] for each φ in `phi`.
    """
    # ensure phi is array
    phi_arr = np.atleast_1d(phi).astype(complex)
    i = 1j

    # Gauss–Legendre nodes & weights on [0,T]
    x, w = leggauss(n_quad)
    s    = 0.5*(x + 1)*T    # shape (n_quad,)
    J    = T/2

    # broadcast s to shape (1, n_quad)
    s_tile = s[None, :]     # (1, n_quad)

    # constants a independent of phi
    a = sigma**2 / 2

    # compute coefficients that depend on phi: shape (M,)
    b = rho * sigma * i*phi_arr - kappa   # (M,)
    c = 0.5 * ((i*phi_arr)**2 - i*phi_arr) # (M,)
    d = np.sqrt(b*b - 4*a*c)               # (M,)
    g = (-b - d)/(-b + d)                  # (M,)

    # expand to broadcast over quadrature nodes: shape (M, n_quad)
    b_tile = b[:, None]
    c_tile = c[:, None]
    d_tile = d[:, None]
    g_tile = g[:, None]

    # closed-form psi2(s) = D(s) at each node: shape (M, n_quad)
    exp_ds = np.exp(-d_tile * s_tile)
    D_s = ((-b_tile - d_tile) / sigma**2) * (1 - exp_ds) / (1 - g_tile * exp_ds)

    # branching function F_s(s) shape (M, n_quad)
    psi1 = i * phi_arr[:, None]
    F_s = (
        0.5*(psi1**2 - psi1)
      + (rho*sigma*psi1 - kappa) * D_s
      + 0.5 * sigma**2 * D_s**2
    )

    xi_flip = xi0(T - s)     
        
    # g0 = v0*np.exp(-kappa*u) + theta*(1 - np.exp(-kappa*u))
    

    # convolution I1 = ∫ F_s(s) * g0(T-s) ds  -> shape (M,)
    I1 = J * np.dot(F_s * xi_flip, w)

    # Why do I need this extra term? W/out it solution does not coincide with classical closed-form ricatti
    # but the math says we do NOT need this term? (TODO)
    I2 = kappa * theta * (J * np.dot(D_s, w))

    # combine and return same shape as input
    exponent = I1 + I2
    phi_out = np.exp(exponent)
    # if input was scalar, return scalar
    return phi_out if phi_out.shape != (1,) else phi_out[0]
```

With

```
xi0      = lambda t: v0*np.exp(-kappa*t) + theta*(1 - np.exp(-kappa*t))
Δ        = 0.01 * np.sqrt(T)          
xi0_bump = lambda t: xi0(t) + (rho * Δ * sigma * np.exp(-kappa * t))
```

Which one can plug into a fast-fourier-transform solver. The resulting implied volatility and strike pairs are then interpolated and evaluated at logmoneyness=$0$. (i.e: $|\beta| = |(iv1 - iv0) / h|$)

I try to replicate Figure 14.3, given below, for the Heston model (purple dots).

The paper details that it is constructed using the following parameters.

\begin{equation} v0=0.117 \\ \kappa=3.37\\ \theta=0.048\\ \sigma=1.99\\ \rho=-0.68\\ \end{equation}

We know that for the Heston, as $T$ goes to $0$, $SSR \approx 2$, and as $T$ goes to $\infty$, $SSR=1$. However, I find the following.

Clearly, as $T$ increases my hypothesis is that the bump becomes too small. i.e $e^{-\kappa u}$ shrinks too fast. If I use just $h=0.01 \times \sqrt{T}$ as the bump, I get the following.

Which is a bit more reaistic, but still not close at all to Figure 14.3

Or am I cutting way too many corners and should I do a proper time discretization? (how?)

And actually, what is interesting, is that I am kinda able to reproduce Figure 7.6

As illustrated below, using Parameter Set II.. :/

So the problems lies somewhere in $\kappa \neq 0 $

## Answer by Rutger Versteegden (score 1)

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

Fixed... (obvious mistake)

```
F_s = (
    0.5*(psi1**2 - psi1)
  + (rho*sigma*psi1 - kappa) * D_s
  + 0.5 * sigma**2 * D_s**2
)
```

Should be

```
F_s = (
    0.5*(psi1**2 - psi1)
  + (rho*sigma*psi1 - 0) * D_s
  + 0.5 * sigma**2 * D_s**2
)
```

And

```
I1+I2
```

Should simply be

```
I1
```

Now I reproduce all figures of Smile Dynamics and Rough Volatility perfectly fine.