Malliavin Monte Carlo Weights for Heston Delta Estimation

# Delta of European Call using Malliavin Calculus


# Delta of European Call using Malliavin Calculus












I am trying to implement the calculation of the Delta under the Heston model using Malliavin Calculus in python. According to "Malliavin Calculus in Finance" by Álos, the Malliavin delta is given as

$$\Delta = \frac{\partial V}{\partial S_0} = \frac{e^{-rT}}{T\sqrt{1-\rho^2}S_0}\mathbb{E}\left[f(S_T)\int_0^T\frac{dB_s}{\sigma_s}\right],$$

where $\rho$ is the correlation between the price- and volatility-process, $T$ is the maturity time of the contract and $B$ is the Brownian motion driving the volatility-process.

The evaluate the integral numerically, I use the following approximation:

$$\int_0^T\frac{dB_s}{\sigma_s}=\sum_{i=0}^N \frac{\sqrt{\Delta t}}{\sigma_t}Z_t,$$

where $N$ is the number of steps, $\Delta t=T/N$ and $Z_t$ is are idd random variables taken from $N(0,1)$ at each time step.

My current attempt for a European call can be found below:

```
# Heston Simumation Parameters for a Black-Scholes Limit
# simulation dependent
S0 = 100.0             # asset price
K = 95                 # strike price
T = 1.0                # time in years
r = 0.02               # risk-free rate
N = 365                # number of time steps in simulation
M = 100000                  # number of simulations
h = 0.01*S0

# Heston dependent parameters
kappa = 0              # rate of mean reversion of variance under risk- 
neutral dynamics
v0 = 0.25**2           # initial variance under risk-neutral dynamics
theta = v0        # long-term mean of variance under risk- 
                           #neutral dynamics
rho = 0.0              # correlation between returns and variances 
                       #under risk-neutral dynamics
xi = 0.0000000001            # volatility of volatility

def heston_sim(S0, K, v0, r, rho, kappa, theta, xi, T, N, M, h):
"""
Inputs:
 - S0, v0: initial parameters for asset and variance
 - rho   : correlation between asset returns and variance
 - kappa : rate of mean reversion in variance process
 - theta : long-term mean of variance process
 - xi : vol of vol / volatility of variance process
 - T     : time of simulation
 - N     : number of time steps
 - M     : number of scenarios / simulations

Outputs:
- asset prices over time (numpy array)
- variance over time (numpy array)
"""

dt = T/N
mu = np.array([0,0])

S = np.full(shape = (N+1, M), fill_value = S0)
S_up = np.full(shape = (N+1, M), fill_value = S0+h)
S_down = np.full(shape = (N+1, M), fill_value = S0-h)
v = np.full(shape = (N+1, M), fill_value = v0)

Z = np.zeros((N+1, M, 2))

#While performing the monte carlo simulation, we want to obtain the malliavin weights simultaneously
delta_weights = np.zeros(M)

for i in range(1,N+1):
    Z1, Z2 = np.random.normal(0, 1, M), np.random.normal(0, 1, M)
    Z_v = Z1
    Z_s = rho * Z1 + np.sqrt(1 - rho**2) * Z2
    Z[i-1, :, 0] = Z_s
    Z[i-1, :, 1] = Z_v
    S[i] = S[i-1] * np.exp((r - 0.5 * v[i-1]) * dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0])
    S_up[i] = S_up[i-1] * np.exp((r - 0.5 * v[i-1]) * dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0])
    S_down[i] = S_down[i-1] * np.exp((r - 0.5 * v[i-1]) * dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0])
    v[i] = np.abs(v[i-1] + kappa * (theta - v[i-1]) * dt + xi * np.sqrt(v[i-1] * dt) * Z[i-1,:,1])
    
    delta_weights += np.array(Z[i-1,:,1]*np.sqrt(dt/v[i-1])) / (T * S0 * np.sqrt(1 - rho**2))

    
ST_up = S_up[-1]
ST_down = S_down[-1]
ST = S[-1]

payout_up = np.maximum(ST_up - K, 0)
payout_down = np.maximum(ST_down - K, 0)
payout = np.maximum(ST - K, 0)

av_payout_up = np.average(payout_up)
av_payout_down = np.average(payout_down)
av_payout = np.average(payout)

#calculate the vanilla greeks via finite difference
delta_findif = (av_payout_up - av_payout_down)/(2*h) * np.exp(-r*T)
gamma_findif = (av_payout_up - 2*av_payout + av_payout_down)/(h**2) * np.exp(-r*T)

delta_mall = np.exp(-r*T) * np.average(np.multiply(np.maximum(ST-K, 0), delta_weights))

return S, v, delta_mall, delta_findif, gamma_findif
```

This code should calculate both the price and volatility process, calculate the finite difference delta/gamma and simultaneously evaluate the integral numerically.

However, when reducing the Heston model to the constant volatility case and comparing the delta with the expected Black-Scholes result, the error is significant, while the Monte Carlo call price gives a very nice match. The theoretical Black-Scholes call price equals 13.43 in comparison to the Monte Carlo price from the code above being 13.402 (this also varies slightly with each time the code is ran). The theoretical delta value equals 0.6591.

I added the finite difference calculation of the delta to the above code which yields a relative stable result of 0.66.

Also, running the code for $M=100,000$ simulations multiple times gives a very different result for the mall delta, while the call price remains stable.

Is there something wrong with my approximation of the integral or my code? Or could it just be that the Malliavin approach for a vanilla call under the Heston model is just highly inaccurate and unstable?

Any help and/or tips are much appreciated!

## Answer by Andrea (score 2, accepted)

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

I am not sure you have copied the right formula, or your notation is not clear.

$$\Delta = \frac{\partial V}{\partial S_0} = \frac{1}{T\sqrt{1-\rho^2}S_0}\mathbb{E}\left[f(S_T)\int_0^T\frac{dB_s}{\sigma_s}\right]$$

What is $dB_s$? You claim: "$B$ is the Brownian motion driving the volatility-process"

I cannot see how this can be true. Compare it to the likelihood ratio (Glassermann p. 403 for instance), and in case of 0 correlation, it will definitely be wrong.

Read this other source:

https://arxiv.org/pdf/2003.01523 (equation 4.53)

It must be the orthogonal part of the Brownian motion driving the Stock, so `Z2` in your notation.

Results are much much better with this change.

Double check the book, and paste the full specification of the SDEs and the Delta weight.

This is my code, which prints as well the MC error for each output. Start with an analysis of how well it performs with simpler payoffs: $1$, then move to $S$, then to an option.

```
import numpy as np
  

def heston_sim(S0, v0, rho, kappa, theta, xi, v_lb, T, N, M, K):
    """
    Inputs:
    - S0, v0: initial parameters for asset and variance
    - rho   : correlation between asset returns and variance
    - kappa : rate of mean reversion in variance process
    - theta : long-term mean of variance process
    - xi : vol of vol / volatility of variance process
    - v_lb : lower bound for the variance
    - T     : time of simulation
    - N     : number of time steps
    - M     : number of scenarios / simulations
    """

    dt = T / N

    # always reuse the same random numbers
    generator = np.random.default_rng(12345)

    S = np.full(shape=(M,), fill_value=S0)
    v = np.full(shape=(M,), fill_value=v0)

    delta_weights = np.full(shape=(M,), fill_value=0.0)

    for i in range(N):  # N number of steps
        Z1 = generator.normal(0, 1, M)
        Z2 = generator.normal(0, 1, M)

        Z_v = Z1
        Z_s = rho * Z1 + np.sqrt(1 - rho**2) * Z2

        S = S * np.exp(-0.5 * v * dt + np.sqrt(v * dt) * Z_s)

        delta_weights += Z2 * np.sqrt(dt / v)

        v = v + kappa * (theta - v) * dt + xi * np.sqrt(v * dt) * Z_v

        # apply a lower bound to avoid "/0"
        # no longer Heston, but maybe more stable
        v = np.maximum(v_lb, v)

    delta_weights /= T * S0 * np.sqrt(1 - rho**2)

    payoff = np.maximum(S - K, 0)

    factor = 1.0 / np.sqrt(M)

    price = np.average(payoff)
    s_price = np.std(payoff) * factor

    # delta
    d_mv_p = np.average(delta_weights * payoff)
    s_mv_p = np.std(delta_weights * payoff) * factor

    # check derivative of "1" (i.e. 0)
    d_mv_1 = np.average(delta_weights * 1)
    s_mv_1 = np.std(delta_weights * 1) * factor

    # check derivative of "s" (i.e. 1)
    d_mv_s = np.average(delta_weights * S)
    s_mv_s = np.std(delta_weights * S) * factor

    return price, s_price, d_mv_p, s_mv_p, d_mv_1, s_mv_1, d_mv_s, s_mv_s

def main():
    S0 = 100
    sigma0 = 0.1
    v0 = sigma0**2
    rho = -0.5
    kappa = 4
    theta = v0
    xi = 0.05
    T = 2.0

    # explicit lower bound in the Heston simulation
    sigma_lb = 0.01  # in vol space
    v_lb = sigma_lb**2

    K = 95
    N = 10  # number of time steps
    M = 1_000_000  # number of paths
    args = S0, v0, rho, kappa, theta, xi, v_lb, T, N, M, K
    out = heston_sim(*args)
    price, s_price, d_mv_p, s_mv_p, d_mv_1, s_mv_1, d_mv_s, s_mv_s  = out
    print(f"{price=}")
    print(f"{s_price=}")
    print(f"{d_mv_p=}")
    print(f"{s_mv_p=}")
    print(f"{d_mv_1=}")  # this should be 0
    print(f"{s_mv_1=}")  # this should be small
    print(f"{d_mv_s=}")  # this should be 1
    print(f"{s_mv_s=}")  # this should be small

    print()

    dS = 0.01  # finite difference delta
    price_down = heston_sim(S0 - dS, v0, rho, kappa, theta, xi, v_lb, T, N, M, K)
    price_up = heston_sim(S0 + dS, v0, rho, kappa, theta, xi, v_lb, T, N, M, K)
    delta_mc = (price_up[0] - price_down[0]) / (2 * dS)
    print(f"{delta_mc=}")  # expected the same as d_mv_p

main()
```