Simulation de mouvement brownien géométrique en Python

Question

J'essaye de simuler le mouvement brownien géométrique en Python, pour fixer le prix d'une option d'appel européen via la simulation Monte-Carlo. Je suis relativement nouveau à Python, et je reçois une réponse que je crois être fausse, car il est loin de converger vers le prix de BS, et les itérations semblent avoir une tendance négative pour une raison quelconque. Toute aide serait appréciée.

import numpy as np 
from matplotlib import pyplot as plt 


S0 = 100 #initial stock price 
K = 100 #strike price 
r = 0.05 #risk-free interest rate 
sigma = 0.50 #volatility in market 
T = 1 #time in years 
N = 100 #number of steps within each simulation 
deltat = T/N #time step 
i = 1000 #number of simulations 
discount_factor = np.exp(-r*T) #discount factor 

S = np.zeros([i,N]) 
t = range(0,N,1) 



for y in range(0,i-1): 
    S[y,0]=S0 
    for x in range(0,N-1): 
     S[y,x+1] = S[y,x]*(np.exp((r-(sigma**2)/2)*deltat + sigma*deltat*np.random.normal(0,1))) 
    plt.plot(t,S[y]) 

plt.title('Simulations %d Steps %d Sigma %.2f r %.2f S0 %.2f' % (i, N, sigma, r, S0)) 
plt.xlabel('Steps') 
plt.ylabel('Stock Price') 
plt.show() 

C = np.zeros((i-1,1), dtype=np.float16) 
for y in range(0,i-1): 
    C[y]=np.maximum(S[y,N-1]-K,0) 

CallPayoffAverage = np.average(C) 
CallPayoff = discount_factor*CallPayoffAverage 
print(CallPayoff) 

Monte-Carlo Simulation Exemple (Stock Prix Simulation)

Je suis actuellement en utilisant Python 3.6.1.

Merci d'avance pour l'aide.

Answer 1

est ici un peu de réécriture de code qui peut rendre la notation de S plus intuitive et vous permettra d'inspecter votre réponse pour le caractère raisonnable.

points initiaux:

• Dans votre code, la deuxième deltat devrait être remplacé par np.sqrt(deltat). Source here (oui, je sais que ce n'est pas le plus officiel, mais les résultats ci-dessous devraient être rassurants).
• Le commentaire concernant la non-annualisation de vos valeurs short et sigma peut être incorrect. Cela n'a rien à voir avec la dérive vers le bas que vous voyez. Vous devez les garder à des taux annualisés. Ceux-ci seront toujours continuellement composés (taux constants).

Tout d'abord, voici une fonction de génération de chemin GBM de Yves Hilpisch - Python Finance, chapter 11. Les paramètres sont expliqués dans le lien mais la configuration est très similaire à la vôtre.

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[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt + 
             sigma * np.sqrt(dt) * rand) 
    return paths 

Définition de vos valeurs initiales (mais en utilisant N=252, le nombre de jours de bourse en 1 an, le nombre d'incréments de temps):

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) 

Générez ensuite les chemins:

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

Maintenant, pour inspecter: paths[-1] vous obtient les valeurs de fin St, à l'expiration:

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

Le paiement, que vous avez maintenant, sera le maximum de (St - K, 0):

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

Si vous tracer ces chemins (facile à utiliser simplement pd.DataFrame(paths).plot(), vous verrez qu'ils » re n'est plus orienté vers le bas, mais que les St sont approximativement log-normalement distribués.

Enfin, voici une vérification raisonnable par 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 

En utilisant des valeurs plus élevées pour i dans la configuration de votre GBM devrait provoquer une plus grande convergence.

Answer 2

Il semble que vous utilisiez la mauvaise formule.

Avez dS_t = S_t (r dt + sigma dW_t) de Wikipedia
Et dW_t ~ Normal(0, dt) de Wikipedia
Alors S_(t+1) = S_t + S_t (r dt + sigma Normal(0, dt))

Je crois donc la ligne devrait plutôt être le suivant:

S[y,x+1] = S[y,x]*(1 + r*deltat + sigma*np.random.normal(0,deltat))