Modelo Heston simulación MC python

Question

Tengo este ejercicio. $\\\\$ Busque valores realistas de los parámetros y calcule el precio de una Call europea con vencimiento $T = 0.5$ y $S_0 = 1$ para los valores de la huelga $K = 0.5,0.6, ......, 1.5$ . Muestra la sonrisa de la volatilidad.

Posteo mi código en python, pero no entiendo por qué no se muestra el smiley de volatilidad.

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from py_vollib_vectorized import vectorized_implied_volatility as implied_vol

# Parameters
# simulation dependent
S0 = 1           # asset price
T = 0.5              # time in years
r = 0.00003            # risk-free rate
N = 125             # number of time steps in simulation
M = 10000         # number of simulations

# Heston dependent parameters
kappa = 1.1              # rate of mean reversion of variance under risk-neutral dynamics
theta = 0.02        # long-term mean of variance under risk-neutral dynamics
v0 = 0.025          # initial variance under risk-neutral dynamics
rho = -0.7              # correlation between returns and variances under risk-neutral dynamics
sigma = 0.1          # volatility of volatility

def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M):
    """
    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
     - sigma : 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)
    """
    # initialise other parameters
    dt = T/N
    mu = np.array([0,0])
    cov = np.array([[1,rho],
                    [rho,1]])

    # arrays for storing prices and variances
    S = np.full(shape=(N+1,M), fill_value=S0)
    v = np.full(shape=(N+1,M), fill_value=v0)

    # sampling correlated brownian motions under risk-neutral measure
    Z = np.random.multivariate_normal(mu, cov, (N,M))

    for i in range(1,N+1):
        S[i] = S[i-1] * np.exp( (r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0] )
        v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0)

    return S, v

rho_p = 0.98
rho_n = -0.98

S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M)
S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M)

fig, (ax1, ax2)  = plt.subplots(1, 2, figsize=(12,5))
time = np.linspace(0,T,N+1)
ax1.plot(time,S_p)
ax1.set_title('Heston Model Asset Prices')
ax1.set_xlabel('Time')
ax1.set_ylabel('Asset Prices')

ax2.plot(time,v_p)
ax2.set_title('Volatility Monte Carlo Simulation')
ax2.set_xlabel('Time')
ax2.set_ylabel('Variance')

plt.show()
# simulate gbm process at time T
gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) )

fig, ax = plt.subplots()

ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax)
ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax)
ax = sns.kdeplot(gbm, label="GBM", ax=ax)

plt.title(r'Asset Price Density under Heston Model')
plt.xlim([0.5, 1.5])
plt.xlabel('$S_T$')
plt.ylabel('Density')
plt.legend()
plt.show()

rho = -0.7
S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M)

# Set strikes and complete MC option price for different strikes
K = np.arange(0.5,1.5)

calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K])

call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore')

plt.plot(K, call_ivs, label=r'Call')

plt.ylabel('Implied Volatility')
plt.xlabel('Strike')

plt.title('Implied Volatility Smile from Heston Model')
plt.legend()
plt.show()

Gracias por el apoyo.

Answer 1

Problemas detectados:

• La cuestión principal: Usted define S0 = 1 que en Python es del tipo Int . La trayectoria de sus acciones se basa en S0 y por lo tanto S será del tipo Int también ( más concretamente, Int32 ). Esto implica que cualquier decimal en sus estimaciones será truncado. Para solucionar esto simplemente definir S0 = 1.0 que ahora es del tipo float .

• Un problema menor: Al definir su conjunto de huelgas $K$ necesitas definir el tamaño del paso (automáticamente fijado en 1) en la función numpy arange, K = np.arange(0.1,1.6,0.1) . Esto le dará el conjunto de huelgas especificado.

• Al probar el código, he detectado que implied_vol función puede se encuentra con algunos problemas de optimización, lo que genera una extraña sonrisa. Conseguí una sonrisa decente (mostrada abajo) después de un par de repeticiones de la función.

Es posible que haya otros pequeños problemas, pero al arreglar los problemas anteriores tengo una solución que funciona.

Verificación:

Me tomé la libertad de arreglar los problemas anteriores y ejecuté el código. Sin embargo, he utilizado los parámetros originales de Heston que se encuentran en el Tutorial en Youtube que ha seguido claramente para producir las simulaciones, la densidad y la sonrisa. Esto da como resultado los siguientes gráficos:

Y la siguiente sonrisa:

Dejaré la interpretación de los gráficos a su criterio. Espero que esto ayude.

---

Apéndice: Código Python

# Parameters
# simulation dependent
S0 = 1.0          # asset price
T = 0.5              # time in years
r = 0.00003            # risk-free rate
#r=0.02
N = 252             # number of time steps in simulation
M = 10000         # number of simulations

# Heston dependent parameters
#kappa = 1.1              # rate of mean reversion of variance under risk-neutral dynamics
#theta = 0.02        # long-term mean of variance under risk-neutral dynamics
#v0 = 0.025          # initial variance under risk-neutral dynamics
#rho = -0.7              # correlation between returns and variances under risk-neutral dynamics
#sigma = 0.1          # volatility of volatility

kappa = 3
theta = 0.20**2
v0=0.25**2
#rho = 0.7
sigma = 0.6 

def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M):
    """
    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
     - sigma : 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)
    """
    # initialise other parameters
    dt = T/N
    mu = np.array([0,0])
    cov = np.array([[1,rho],
                    [rho,1]])

    # arrays for storing prices and variances
    S = np.full(shape=(N+1,M), fill_value=S0)
    v = np.full(shape=(N+1,M), fill_value=v0)

    # sampling correlated brownian motions under risk-neutral measure
    Z = np.random.multivariate_normal(mu, cov, (N,M))

    for i in range(1,N+1):
        S[i] = S[i-1] *  np.exp((r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0])
        v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0)

    return S, v

rho_p = 0.98
rho_n = -0.98

S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M)
S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M)

fig, (ax1, ax2)  = plt.subplots(1, 2, figsize=(12,5))
time = np.linspace(0,T,N+1)
ax1.plot(time,S_p)
ax1.set_title('Heston Model Asset Prices')
ax1.set_xlabel('Time')
ax1.set_ylabel('Asset Prices')

ax2.plot(time,v_p)
ax2.set_title('Volatility Monte Carlo Simulation')
ax2.set_xlabel('Time')
ax2.set_ylabel('Variance')

plt.show()
# simulate gbm process at time T
gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) )

fig, ax = plt.subplots()

ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax)
ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax)
ax = sns.kdeplot(gbm, label="GBM", ax=ax)

plt.title(r'Asset Price Density under Heston Model')
plt.xlim([0.5, 1.5])
plt.xlabel('$S_T$')
plt.ylabel('Density')
plt.legend()
plt.show()

rho = -0.7
S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M)

# Set strikes and complete MC option price for different strikes

K = np.arange(0.1,1.6,0.1)

calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K])

call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore')

plt.plot(K, call_ivs, label=r'Call')

plt.ylabel('Implied Volatility')
plt.xlabel('Strike')

plt.title('Implied Volatility Smile from Heston Model')
plt.legend()
plt.show()