Simulación numérica del modelo de Heston

Question

Estoy tratando de simular en Python trayectorias aleatorias para el precio de un activo general como lo describe el modelo de Heston:

\begin{equation} \begin{aligned} dS_t &= \mu S_t dt + \sqrt{\nu_t} S_t dW^S_t \\ d\nu_t &= \kappa(\theta - \nu_t) dt + \xi \sqrt{\nu_t} dW^{\nu}_t \\ \textrm{Corr}[W^S_t, W^{\nu}_t] &= \rho \end{aligned} \end{equation}

donde:

• $W_t^{S}$ y $W_t^{\nu}$ son dos movimientos brownianos estándar con una correlación $\rho$.
• $\nu _{t}$ es la varianza instantánea.
• $\mu$ es la tasa de retorno del activo.
• $\theta$ es la varianza a largo plazo.
• $\kappa$ es la tasa a la que $\nu_t$ vuelve a $\theta$.
• $\xi$ es la volatilidad de la volatilidad instantánea.

Por lo tanto, implementé la siguiente función:

def HeMC (S0, mu, v0, rho, kappa, theta, xi, T, dt):

    # Generar una simulación de Monte Carlo para el modelo de Heston

    # Generar un Movimiento Browniano aleatorio
    MU  = np.array([0, 0])
    COV = np.matrix([[1, rho], [rho, 1]])
    W   = np.random.multivariate_normal(MU, COV, T)
    W_S = W[:,0]
    W_v = W[:,1]

    # Generar trayectorias
    vt    = np.zeros(T)
    vt[0] = v0
    St    = np.zeros(T)
    St[0] = S0
    for t in range(1,T):
        vt[t] = np.abs(vt[t-1] + kappa*(theta-np.abs(vt[t-1]))*dt + xi*np.sqrt(np.abs(vt[t-1]))*W_v[t])
        St[t] = St[t-1]*np.exp((mu - 0.5*vt[t])*dt + np.sqrt(vt[t]*dt)*W_S[t])

    return St, vt

El problema es que cuando ejecuto esta función con los siguientes parámetros, a menudo obtengo trayectorias que parecen no tener sentido. Especialmente la volatilidad instantánea no parece revertir a la media, sino que a menudo sigue trayectorias salvajes.

T     = 252
dt    = 1/252
S0    = 100 # Precio inicial
mu    = 0.1 # Retorno esperado
sigma = 0.2 # Volatilidad
rho   = -0.2 # Correlación
kappa = 0.3 # Tasa de reversión
theta = 0.2 # Volatilidad a largo plazo
xi    = 0.2 # Volatilidad de la volatilidad instantánea
v0    = 0.2 # Volatilidad instantánea inicial

Creo que el problema está en la forma en que discreticé el proceso o tal vez hay algunos errores en mi código, pero no pude encontrar nada.

Gracias por tu ayuda.

Answer 1

Hay dos errores en el código:

1) En la línea

vt[t] = np.abs(vt[t-1] + kappa*(theta-np.abs(vt[t-1]))*dt + xi*np.sqrt(np.abs(vt[t-1]))*W_v[t]*np.sqrt(dt])

olvidaste multiplicar W_v[t] por np.sqrt(dt).
Esta es la razón por la que la volatilidad aumenta tanto.

2) La línea

St[t] = St[t-1]*np.exp((mu - 0.5*vt[t])*dt + np.sqrt(vt[t]*dt)*W_S[t])

debería ser

St[t] = St[t-1]*np.exp((mu - 0.5*vt[t-1])*dt + np.sqrt(vt[t-1]*dt)*W_S[t])

Tampoco es necesario usar tres veces la función np.abs. Una es suficiente. (La más externa).

Answer 2

El proceso de varianza en el modelo Heston (es decir, proceso CIR) es notorio para la simulación. Un esquema típico de Euler/Milstein termina con varianza negativa en una gran parte de las trayectorias. Hacerla positiva con np.abs() distorsiona la verdadera distribución.

Considere el esquema QE de Andersen (2008) para resolver este problema:

• Andersen L (2008) Simulación simple y eficiente del modelo de volatilidad estocástica de Heston. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189

Alternativamente, considere la distribución ji-cuadrada no central para simular vt, que es la propiedad analítica del proceso CIR. El código a continuación es solo una idea. (Por favor modifíquelo adecuadamente si hay errores de sintaxis.)

chi_df = 4 * theta * kappa / xi**2
exph = np.exp(-kappa*dt/2)
phi = 4*kappa / xi**2 / (1/exph - exph)
chi_nonc = vt[t-1] * exph * phi
vt[t] = (exph / phi) * np.random.noncentral_chisquare(df=chi_df, nonc=chi_nonc)