2 votos

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

Estoy tratando de realizar Simulaciones de Monte Carlo usando números normales estándar casi 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 poco realistas de la trayectoria de los 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):
    """
    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(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 ser generados, 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

Cualquier ayuda es apreciada.

0 votos

¿Qué es una aplicación "sencilla"? ¿Qué has encontrado y por qué no está bien? ¿Implementación en pseudocódigo o en qué lenguaje?

0 votos

@gg: He editado el post para incluir el funcionamiento

1 votos

¿Así que tal vez "sobol_seq.i4_sobol_generate" es incorrecto? ¿Qué has hecho para comprobar la corrección? ¿Y cómo podemos ver en tu código que "8.11596295e-02 no debe ser generado"?

3voto

Esto sucede porque estás utilizando los mismos números aleatorios (psuedo/cuasi) para cada paso de tiempo.

en su 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) se están configurando los números aleatorios para que sean los mismos en cada punto de tiempo. El impacto de esto es que un camino que comienza con una probabilidad improbable obtendrá esa misma probabilidad improbable en cada paso de tiempo, y así todos sus caminos divergen:

enter image description here

Tienes que 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 su código, ver 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 se ven así:

enter image description here

o con más caminos y alguna transparencia añadida: enter image description here

Ahora, hay un problema con esto -> sobol_seq se limita a 40 dimensiones. Hay dos maneras de arreglar esto:

  1. Escriba su propio código de secuencia sobol e incluya más direction numbers para poder ir a dimensiones superiores. Esto no es trivial, pero de ninguna manera se puede lograr. Existe un código aquí, crédito a Stephen Joe, y Frances Kuo con números de dirección para llevarle hasta las dimensiones de 21k.
  2. Se utiliza un Puente Browniano . Utilizando esta técnica, se puede reatribuir la naturaleza de baja discrepancia de la secuencia de sobol a los puntos que controlan la mayor parte de la varianza en el MC. A los efectos de generar simplemente trayectorias, esto significa que se utiliza la primera dimensión para el paso de tiempo total, luego el segundo punto para generar el punto puente en el medio, luego se hace recursivamente lo mismo y se rellenan los huecos. una vez que se agotan los números sobol, se ha dividido el espacio en 40 subperíodos de tiempo. A continuación, simplemente rellena los huecos de estos espacios con números aleatorios, ya que las ganancias que obtendrá ahora por la baja discrepancia son extremadamente pequeñas.

0 votos

Para la generación de secuencias Sobol en altas dimensiones, un colega ha hecho un wrapper en C y Python para la generación con hace uso del Intel MKL, que él mantiene: bitbucket.org/croci/mkl_sobol/src/master

0 votos

@will Gracias por tu respuesta y detalles. ¿Puede usted por favor ampliar en el puente browniano para llenar los números entre los de Sobol? si todo lo que necesito es uniforme números aleatorios de baja discrepancia con una dimensión muy alta (digamos 500), cómo puedo generar números a partir de la secuencia de sobol (digamos 10 dimensiones) y luego llenar el resto de las dimensiones. ¿Puede usted por favor enumerar los pasos y también hacer referencia a cualquier enlace? ¿Existe un puente browniano para variables aleatorias no formadas o sólo para variables aleatorias normales?

0 votos

@toing ¿La varianza aportada por cada uno de tus números aleatorios es la misma? o ¿puedes asignar una cantidad diferente de varianza a cada variable requerida? Si quieres más, entonces puedes mirar este código de S. Joe y F.Y.Kuo, tiene números de dirección para números de sobol hasta algo más de 21k.

1voto

oliversm Puntos 515

Si utiliza puntos de baja discrepancia, debería aleatorizarlos antes de transformarlos en una distribución normal. Hay muchas maneras de hacer esto, pero para simplificar la codificación he utilizado una traducción uniforme ( % 1 en Python). Como ejemplo, esto se vería como

import numpy as np
from scipy.stats import norm
import sobol_seq as ss
randomisation = np.random.random()  # Ensures the sequence is randomised.
uniforms = np.concatenate(ss.i4_sobol_generate(1, 30))  # Some definitions start with 0!
uniforms = (randomisation + uniforms) % 1  # Digital shifting would also be an option
normals = norm.ppf(uniforms)  # The inverse transform method.

En cuanto a los valores pequeños, dado que el proceso logarítmico de un movimiento browniano geométrico $S_t$ es un movimiento browniano a la deriva con media $\log(S_0) + (\mu - \tfrac{\sigma^2}{2})t$ y la varianza $\sigma^2 t$ puedes calcular la probabilidad de que el mínimo de este proceso caiga por debajo de un nivel determinado. En tu caso unos 4 órdenes de magnitud en base-10 y unos 9 en base- $\rm{e}$ ). La expresión de 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.

0 votos

Por favor, ¿puede darme una referencia que explique por qué es necesario aleatorizar los puntos antes de utilizarlos? No veo la necesidad de hacerlo.

1 votos

@will Esto se conoce como aleatorio cuasi-Monte Carlo . Como la secuencia es determinista y la estimación de la primera mitad de la secuencia no es independiente de la segunda mitad, no podemos utilizar la CLT. En su lugar, lo único que podemos utilizar para limitar el error es la desigualdad de Koksma-Hlawka, que es difícil de estimar. Para evitarlo, utilizamos estimadores aleatorios. Se puede encontrar buen material sobre esto en "Randomized Quasi-Monte Carlo: An Introduction for Practitioners" (L'Ecuyer) y "Quasi-Monte Carlo Sampling" (Owen). Si quiere recuperar un intervalo de confianza, la aleatorización lo conseguirá.

0voto

tjhorner Puntos 8

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

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

Finanhelp.com

FinanHelp es una comunidad para personas con conocimientos de economía y finanzas, o quiere aprender. Puedes hacer tus propias preguntas o resolver las de los demás.

Powered by:

X