As before, translation into Python and NumPy code is straightforward:
In [ 17 ]: I = 10000
M = 50
dt = T / M
S = np.zeros((M + 1 , I))
S[ 0 ] = S0
for t in range( 1 , M + 1 ):
S[t] = S[t - 1 ] * np.exp((r - 0.5 * sigma ** 2 ) * dt
+ sigma * np.sqrt(dt) * npr.standard_normal(I))
The resulting end values for the index level are log-normally distributed again, as
Figure 10-5 illustrates:
In [ 18 ]: plt.hist(S[- 1 ], bins= 50 )
plt.xlabel(‘index level’)
plt.ylabel(‘frequency’)
plt.grid(True)
Figure 10-5. Simulated geometric Brownian motion at maturity
The first four moments are also quite close to those resulting from the static simulation
approach:
In [ 19 ]: print_statistics(S[- 1 ], ST2)
Out[19]: statistic data set 1 data set 2
–––––––––––––––
size 10000.000 10000.000
min 25.531 27.266
max 425.051 358.997
mean 110.900 110.528
std 40.135 40.894
skew 1.086 1.150
kurtosis 2.224 2.273
Figure 10-6 shows the first 10 simulated paths:
In [ 20 ]: plt.plot(S[:, : 10 ], lw=1.5)
plt.xlabel(‘time’)
plt.ylabel(‘index level’)
plt.grid(True)