For what follows we need a number of Python libraries, including scipy.stats and
statsmodels.api:
In [ 1 ]: import numpy as np
np.random.seed( 1000 )
import scipy.stats as scs
import statsmodels.api as sm
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
Let us define a function to generate Monte Carlo paths for the geometric Brownian motion
(see also Chapter 10):
In [ 2 ]: def gen_paths(S0, r, sigma, T, M, I):
”’ Generates Monte 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 = 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)
rand = (rand - rand.mean()) / rand.std()
paths[t] = paths[t - 1 ] * np.exp((r - 0.5 * sigma ** 2 ) * dt +
sigma * np.sqrt(dt) * rand)
return paths
The following is a possible parameterization for the Monte Carlo simulation, generating,
in combination with the function gen_paths, 250,000 paths with 50 time steps each:
In [ 3 ]: S0 = 100.
r = 0.05
sigma = 0.2
T = 1.0
M = 50
I = 250000
In [ 4 ]: paths = gen_paths(S0, r, sigma, T, M, I)
Figure 11-1 shows the first 10 simulated paths from the simulation:
In [ 5 ]: plt.plot(paths[:, : 10 ])
plt.grid(True)
plt.xlabel(‘time steps’)
plt.ylabel(‘index level’)
Our main interest is in the distribution of the log returns. The following code generates an
ndarray object with all log returns:
In [ 6 ]: log_returns = np.log(paths[ 1 :] / paths[ 0 :- 1 ])
