幾何ブラウン運動のシミュレーション

dS = aSdt + bdW

の確率微分方程式の解は以下になる。

S = S0 exp((r - 1/2 σ**)t + σW)

この式で注意するポイントは３つある。

1点めは基本は金利rで増えること

2点目はボラティリティが大きいと金利rで増えた分が減少すること

3点目は不確定のノイズWがボラティリティと掛け合わせで影響を与えること

実行結果

Pythonコード

import math
import numpy as np
import scipy.stats as scs
import statsmodels.api as sm
from pylab import mpl, plt

plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
%matplotlib inline

def gen_paths(S0, r, sigma, T, M, I):
   ''' Generate dMonte Carlo paths for geometric Brownian motion.
   
   Parameters
   ==========
   S0: float
       initial stock/index value
   r: float
       constant short rate
   sigma: float
       constant volatility
   T: float
       final time horizon
   M: int
       number of time steps/intervals
   I: int
       number of paths to be simulated
       
   Returns
   =======
   paths: ndarray, shape (M + 1, I)
       simulated paths given the parameters
   '''
   dt = T / M
   paths = np.zeros((M + 1, I))
   paths[0] = S0
   for t in range(1, M + 1):
       rand = np.random.standard_normal(I)
       rand = (rand - rand.mean()) / rand.std()  
       paths[t] = paths[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt +
                                        sigma * math.sqrt(dt) * rand)  
   return paths
   S0 = 100.  
r = 0.05  
sigma = 0.3  
T = 1.0   
M = 250  
I = 100 
np.random.seed(1000)

paths = gen_paths(S0, r, sigma, T, M, I)

plt.figure(figsize=(10, 6))
plt.plot(paths[:, :10000])
plt.xlabel('time steps')
plt.ylabel('index level');