Simulación de movimiento browniano geométrico en Python

Solución:

Aquí hay un poco de reescritura de código que puede hacer la notación de S más intuitivo y le permitirá inspeccionar su respuesta para ver si es razonable.

Puntos iniciales:

• En tu código, el segundo deltat debe ser reemplazado por np.sqrt(deltat). Fuente aquí (sí, sé que no es la más oficial, pero los resultados a continuación deberían ser tranquilizadores).
• El comentario sobre la anulación de la anualización de sus valores de tasa corta y sigma puede ser incorrecto. Esto no tiene nada que ver con la deriva descendente que está viendo. Debe mantenerlos a tasas anualizadas. Estos siempre serán tasas de capitalización continua (constantes).

Primero, aquí hay una función de generación de ruta de GBM de Yves Hilpisch: Python para finanzas, capítulo 11. Los parámetros se explican en el enlace pero la configuración es muy similar a la suya.

def gen_paths(S0, r, sigma, T, M, I):
    dt = float(T) / M
    paths = np.zeros((M + 1, I), np.float64)
    paths[0] = S0
    for t in range(1, M + 1):
        rand = np.random.standard_normal(I)
        paths
                                         sigma * np.sqrt(dt) * rand)
    return paths

Estableciendo sus valores iniciales (pero usando N=252, número de días de negociación en 1 año, como el número de incrementos de tiempo):

S0 = 100.
K = 100.
r = 0.05
sigma = 0.50
T = 1
N = 252
deltat = T / N
i = 1000
discount_factor = np.exp(-r * T)

Luego genera las rutas:

np.random.seed(123)
paths = gen_paths(S0, r, sigma, T, N, i)

Ahora, para inspeccionar: paths[-1] te da el final St valores, al vencimiento:

np.average(paths[-1])
Out[44]: 104.47389541107971

La recompensa, como la tiene ahora, será el máximo de (St - K, 0):

CallPayoffAverage = np.average(np.maximum(0, paths[-1] - K))
CallPayoff = discount_factor * CallPayoffAverage
print(CallPayoff)
20.9973601515

Si traza estos caminos (fácil de usar pd.DataFrame(paths).plot(), verá que ya no tienen una tendencia a la baja, pero que el Sts tienen una distribución aproximadamente logarítmica normal.

Por último, aquí hay una verificación de cordura a través de BSM:

class Option(object):
    """Compute European option value, greeks, and implied volatility.

    Parameters
    ==========
    S0 : int or float
        initial asset value
    K : int or float
        strike
    T : int or float
        time to expiration as a fraction of one year
    r : int or float
        continuously compounded risk free rate, annualized
    sigma : int or float
        continuously compounded standard deviation of returns
    kind : str, {'call', 'put'}, default 'call'
        type of option

    Resources
    =========
    http://www.thomasho.com/mainpages/?download=&act=model&file=256
    """

    def __init__(self, S0, K, T, r, sigma, kind='call'):
        if kind.istitle():
            kind = kind.lower()
        if kind not in ['call', 'put']:
            raise ValueError('Option type must be 'call' or 'put'')

        self.kind = kind
        self.S0 = S0
        self.K = K
        self.T = T
        self.r = r
        self.sigma = sigma

        self.d1 = ((np.log(self.S0 / self.K)
                + (self.r + 0.5 * self.sigma ** 2) * self.T)
                / (self.sigma * np.sqrt(self.T)))
        self.d2 = ((np.log(self.S0 / self.K)
                + (self.r - 0.5 * self.sigma ** 2) * self.T)
                / (self.sigma * np.sqrt(self.T)))

        # Several greeks use negated terms dependent on option type
        # For example, delta of call is N(d1) and delta put is N(d1) - 1
        self.sub = {'call' : [0, 1, -1], 'put' : [-1, -1, 1]}

    def value(self):
        """Compute option value."""
        return (self.sub[self.kind][1] * self.S0
               * norm.cdf(self.sub[self.kind][1] * self.d1, 0.0, 1.0)
               + self.sub[self.kind][2] * self.K * np.exp(-self.r * self.T)
               * norm.cdf(self.sub[self.kind][1] * self.d2, 0.0, 1.0))
option.value()
Out[58]: 21.792604212866848

Usando valores más altos para i en su configuración de GBM debería causar una convergencia más cercana.

---

---

Utiliza Nuestro Buscador