Las simulaciones de Monte Carlo en Python utilizando números normales estándar cuasi-aleatorios usando secuencias sobol dan valores erróneos

Question

Estoy tratando de realizar simulaciones de Monte Carlo utilizando números normales estándar cuasi aleatorios. Entiendo que podemos usar secuencias de Sobol para generar números uniformes, y luego usar la transformación integral de probabilidad para convertirlos en números normales estándar. Mi código da valores no realistas de la trayectoria de activos simulados:

import sobol_seq
import numpy as np
from scipy.stats import norm

def i4_sobol_generate_std_normal(dim_num, n, skip=1):
    """
    Genera variables cuasi-aleatorias multivariadas estándar normales.
    Parámetros:
      Entrada, entero dim_num, la dimensión espacial.
      Entrada, entero n, el número de puntos a generar.
      Entrada, entero SKIP, el número de puntos iniciales a omitir.
      Salida, np array real de forma (n, dim_num).
    """

    sobols = sobol_seq.i4_sobol_generate(dim_num, n, skip)

    normals = norm.ppf(sobols)

    return normals

def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):
    dt = float(Ttm) / TradingDaysInAYear
    paths = np.zeros((TradingDaysInAYear + 1, NoOfPaths), np.float64)
    paths[0] = UnderlyingPrice
    for t in range(1, TradingDaysInAYear + 1):
        rand = i4_sobol_generate_std_normal(1, NoOfPaths)
        lRand = []
        for i in range(len(rand)):
            a = rand[i][0]
            lRand.append(a)
        rand = np.array(lRand)

        paths[t] = paths[t - 1] * np.exp((RiskFreeRate - 0.5 * Volatility ** 2) * dt + Volatility * np.sqrt(dt) * rand)
    return paths

GBM(1, 252, 8, 100., 0.05, 0.5)

array([[1.00000000e+02, 1.00000000e+02, 1.00000000e+02, ...,
        1.00000000e+02, 1.00000000e+02, 1.00000000e+02],
       [9.99702425e+01, 1.02116774e+02, 9.78688323e+01, ...,
        1.00978615e+02, 9.64128959e+01, 9.72154915e+01],
       [9.99404939e+01, 1.04278354e+02, 9.57830834e+01, ...,
        1.01966807e+02, 9.29544649e+01, 9.45085180e+01],
       ...,
       [9.28295879e+01, 1.88049044e+04, 4.58249200e-01, ...,
        1.14117599e+03, 1.08089096e-02, 8.58754653e-02],
       [9.28019642e+01, 1.92029616e+04, 4.48483141e-01, ...,
        1.15234371e+03, 1.04211828e-02, 8.34842557e-02],
       [9.27743486e+01, 1.96094448e+04, 4.38925214e-01, ...,
        1.16362072e+03, 1.00473641e-02, 8.11596295e-02]])

Valores como 8.11596295e-02 no deberían generarse, por lo que creo que hay algo mal en el código.

Referencias: https://stats.stackexchange.com/questions/27450/best-method-for-transforming-low-discrepancy-sequence-into-normal-distribution, https://stackoverflow.com/questions/9412339/recommendations-for-low-discrepancy-e-g-sobol-quasi-random-sequences-in-pytho, https://github.com/naught101/sobol_seq

Se agradece cualquier ayuda.

Answer 1

Está ocurriendo porque estás utilizando los mismos números aleatorios (pseudo/cuasi) para cada paso de tiempo.

en tu código aquí:

def GBM(Ttm, TradingDaysInAYear, NoOfPaths, UnderlyingPrice, RiskFreeRate, Volatility):
    dt = float(Ttm) / TradingDaysInAYear
    paths = np.zeros((TradingDaysInAYear + 1, NoOfPaths), np.float64)
    paths[0] = UnderlyingPrice
    for t in range(1, TradingDaysInAYear + 1):
        rand = i4_sobol_generate_std_normal(1, NoOfPaths) # <-- this is the issue
        lRand = []
        for i in range(len(rand)):
            a = rand[i][0]
            lRand.append(a)
        rand = np.array(lRand)

        paths[t] = paths[t - 1] * np.exp((RiskFreeRate - 0.5 * Volatility ** 2) * dt + Volatility * np.sqrt(dt) * rand)
    return paths

en esta línea rand = i4_sobol_generate_std_normal(1, NoOfPaths) estás estableciendo que los números aleatorios sean iguales en cada punto de tiempo. El impacto de esto es que una ruta que comienza con una probabilidad improbable obtendrá esa misma probabilidad improbable en cada paso de tiempo, y por lo tanto todas tus rutas divergen:

Necesitas generar todos los números aleatorios que necesitas para todas tus diferentes dimensiones, donde cada paso de tiempo cuenta como una dimensión diferente. He reescrito tu código, mira aquí:

import sobol_seq
import numpy as np
from scipy.stats import norm

from matplotlib import pyplot

def i4_sobol_generate_std_normal(dim_num, n, skip=1):
    """
    Generates multivariate standard normal quasi-random variables.
    Parameters:
      Input, integer dim_num, the spatial dimension.
      Input, integer n, the number of points to generate.
      Input, integer SKIP, the number of initial points to skip.
      Output, real np array of shape (n, dim_num).
    """
    sobols = sobol_seq.i4_sobol_generate(dim_num, n, skip)
    normals = norm.ppf(sobols)
    return normals

def GBM(s0, r, vol, t, n_paths=8192, dt_days=3, sbol_skip=0):
    dt = dt_days / 365.25

    time_points = []
    t_ = 0
    while t_ + dt < t:
        t_ += dt
        time_points.append(t_)

    time_points.append(t)
    time_points = np.array(time_points)

    n_time_points = len(time_points)

    rand = i4_sobol_generate_std_normal(n_time_points, n_paths, skip=sbol_skip)

    paths = np.zeros((n_paths, n_time_points))
    paths[:,0] = s0
    for t_i in range(1, n_time_points):
        dt_ = time_points[t_i] - time_points[t_i-1]
        paths[:, t_i] = paths[:, t_i-1] * np.exp((r - 0.5*vol**2) * dt_ + np.sqrt(dt_) * vol * rand[:, t_i])

    return time_points, paths

n_paths = 1024
time_points, paths = GBM(100, 0, 0.2, 0.3, n_paths=n_paths)

fig = pyplot.figure()
ax = fig.add_subplot(1,1,1)

for i in range(n_paths):
    ax.plot(time_points, paths[i, :], c=(1,0,0,0.05))

ax.legend()

pyplot.show()

y añadí algo de código para trazar las rutas, ahora lucen así:

o con más rutas y algo de transparencia añadida:

Ahora, hay un problema con esto -> sobol_seq está limitado a 40 dimensiones. Hay 2 formas de solucionar esto:

1. Escribe tu propio código de secuencia de sobol e incluye más números de dirección para que puedas ir a dimensiones más altas. Esto no es trivial, pero de ninguna manera inalcanzable. Hay un código aquí, créditos a Stephen Joe y Frances Kuo, completo con números de dirección que te llevarán hasta 21k dimensiones.
2. Utiliza un Puente Browniano. Usando esta técnica, puedes redistribuir la naturaleza de baja discrepancia de la secuencia de sobol a los puntos que controlan la mayor varianza en el MC. Para los propósitos de simplemente generar rutas, esto significa que utilizas la primera dimensión para el paso de tiempo total, luego el segundo punto para generar el punto de puente en el medio, luego haces lo mismo de forma recursiva y llenas los espacios en blanco. una vez que te quedas sin números de sobol, has dividido el espacio en 40 subperíodos. Luego simplemente llenas los espacios en blanco de estos espacios con números aleatorios, ya que las ganancias que obtendrás ahora por la baja discrepancia son extremadamente pequeñas.

Answer 2

Como sugerencia, si estás utilizando puntos de discrepancia baja, entonces deberías aleatorizarlos antes de transformarlos a la distribución normal. Hay muchas maneras de hacer esto, pero para simplificar el código he utilizado una traslación uniforme (% 1 en Python). Como ejemplo, esto se vería así

import numpy as np
from scipy.stats import norm
import sobol_seq as ss
randomisation = np.random.random()  # Asegura que la secuencia esté aleatorizada.
uniforms = np.concatenate(ss.i4_sobol_generate(1, 30))  # ¡Algunas definiciones comienzan con 0!
uniforms = (randomisation + uniforms) % 1  # La traslación digital también sería una opción
normals = norm.ppf(uniforms)  # El método de transformación inversa.

En cuanto a los valores pequeños, dado que el log-proceso de un movimiento Browniano geométrico $S_t$ es un movimiento Browniano con deriva con media $\log(S_0) + (\mu - \tfrac{\sigma^2}{2})t$ y varianza $\sigma^2 t$ puedes calcular qué tan probable es que el mínimo de este proceso caiga por debajo de un cierto nivel. En tu caso alrededor de 4 órdenes de magnitud en base 10 y alrededor de 9 en base-$\rm{e}$). La expresión para esta probabilidad se puede encontrar aquí. Puedes evaluar esto para los parámetros que tienes y ver con qué frecuencia esperas ver un valor tan pequeño, y compararlo con tus simulaciones.

Answer 3

Ahora puedes usar pytorch que puede generar números de sobol mezclados (números aleatorios de baja discrepancia) en un QMC.

https://pytorch.org/docs/stable/generated/torch.quasirandom.SobolEngine.html