Geometrische Brownsche Bewegungssimulation in Python

Question

Ich versuche, Geometrische Brownsche Bewegung in Python zu simulieren, um eine Europäische Call-Option durch Monte-Carlo-Simulation zu berechnen. Ich bin relativ neu in Python, und ich bekomme eine Antwort, die ich glaube falsch zu sein, da es bei weitem nicht zum BS-Preis konvergiert, und die Iterationen scheinen aus irgendeinem Grund negativ zu sein. Jede Hilfe wäre willkommen.

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 Beispiel (Aktienkurs Simulation)

ich zur Zeit bin mit Python 3.6.1.

Vielen Dank im Voraus für die Hilfe.

Answer 1

Hier ist ein bisschen Umschreiben des Codes, der kann Machen Sie die Notation S intuitiver und können Sie Ihre Antwort auf Angemessenheit prüfen.

Anfangspunkte:

• In Ihrem Code, die zweiten deltat sollen durch np.sqrt(deltat) ersetzt werden. Quelle here (ja, ich weiß, es ist nicht die offizielle, aber die Ergebnisse unten sollten beruhigend sein).
• Der Kommentar bezüglich der Annularisierung Ihrer Werte für kurze Rate und Sigma ist möglicherweise falsch. Das hat nichts mit der Abwärtsdrift zu tun, die du siehst. Sie müssen diese zu annualisierten Preisen behalten. Diese werden immer kontinuierlich (konstante) Raten sein.

Als erstes ist hier ein GBM-Pfad-Funktion von Yves Hilpisch zu erzeugen - Python für Finanzen, chapter 11. Die Parameter sind im Link erklärt, aber das Setup ist sehr ähnlich zu Ihrem.

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 

Einstellen Anfangswerte (aber mit N=252, Anzahl der Handelstage in 1 Jahr, wie die Anzahl der Zeitschritte):

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) 

dann die Pfade erzeugen:

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

Jetzt zu inspizieren: paths[-1] bekommt Sie die Endung St Werte, bei Ablauf:

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

Die Auszahlung, wie Sie jetzt haben, wird das Maximum von (St - K, 0) sein:

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

Wenn Sie diese Pfade (leicht plotten nur pd.DataFrame(paths).plot() verwenden, werden Sie sehen, dass sie‘ Sie sind nicht mehr abwärts gerichtet, aber die St sind ungefähr logarithmisch normalverteilt.

Schließlich ist hier eine Plausibilitätsprüfung durch 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 

höhere Werte für i in Ihrem GBM Setup näher Konvergenz führen sollte.

Answer 2

Es sieht so aus, als ob Sie die falsche Formel verwenden.

Haben dS_t = S_t (r dt + sigma dW_t) von Wikipedia
Und dW_t ~ Normal(0, dt) von Wikipedia
So S_(t+1) = S_t + S_t (r dt + sigma Normal(0, dt))

So glaube ich, die Linie anstatt dies sein sollte:

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