I have run some stock simulation using GBM. I have written two functions, one is using the for-loop and is more computationally expensive. After that, I wrote another which is vectorized and thus more efficient. However, after displaying the results, it seems that the functions display a bit different behaviour.Even though the plots look very similar the mean is different for both of them - the first function has a mean around 105 and the other around 102. I do not think that such a difference is just due to randomness (especially when run with 1,000 simulations)
import numpy as np
import matplotlib.pyplot as plt
import time
np.random.seed(123)
def MC_gbm1(S0, mu, sigma, T, M, n):
start_time2 = time.time()
dt = float(T) / n
paths = np.zeros((n+1,M), np.float64)
paths[0] = S0
for t in range(1, n + 1):
rand = np.random.normal(0,np.sqrt(dt), M)
paths[t] = paths[t-1] * np.exp((mu- 0.5 * sigma ** 2) * dt +
sigma * rand)
comp_time2 = round(time.time() - start_time2,4)
print("Computation time:", comp_time2)
return paths
def MC_gbm2(S0, mu, sigma, T, M, n):
#timesteps = np.linspace(0, T, n+1)[:,None]
dt = T / n
dW = np.random.normal(0, np.sqrt(dt), size=(n, M))
dW_cumsum = np.cumsum(dW, axis=0)
paths = S0 * np.exp((mu - 0.5 * sigma**2) * dt + sigma * dW_cumsum)
result = np.concatenate((np.full((1,M), S0), paths), axis=0)
return result
# Parameters
mu = 0.1
n = 100 # number of steps
T = 1 # time in years
M = 1000 #no of sims
S0 = 100
sigma = 0.3
# GBM Simulation 1
St = MC_gbm1(S0, mu, sigma, T, M, n)
# GBM Simulation 2
St3 = MC_gbm2(S0, mu, sigma, T, M, n)
# Plot the results
sequence = np.linspace(0, T, n+1)
tt = np.full(shape = (M,n+1), fill_value=sequence).T
plt.plot(tt,St)
plt.title("GBM simulation 1")
plt.ylabel("Stock price")
plt.xlabel("Years")
plt.show()
plt.plot(tt,St2)
plt.title("GBM simulation 2")
plt.ylabel("Stock price")
plt.xlabel("Years")
plt.show()
enter image description here enter image description here
I have tried to modify the functions in several ways but nothing really changed. Ideally, I would like to see two very means and plots (of course taking into account the slight randomness given by the stochastic nature)
In order to understand the difference, note that resetting the seed between the two following operations results in the exact same matrix of noise:
n = 10
M = 1000
dt = 1 / 10
rands = np.empty((n, M))
# Operation 1
np.random.seed(123)
for i in range(n):
rands[i] = np.random.normal(0, np.sqrt(dt), M)
# Operation 2
np.random.seed(123)
dW = np.random.normal(0, np.sqrt(dt), size=(n, M))
print("Equal matrices?", np.allclose(rands, dW)) # True
So the difference must come from the order and amount of operations in each function. In fact, if you turn off the random noise in each function by doing
paths[t] = paths[t-1] * np.exp((mu- 0.5 * sigma ** 2) * dt + sigma * np.zeros_like(rand))
and
paths = S0 * np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.zeros_like(dW_cumsum))
both functions produce different matrices:
# GBM Simulation 1
np.random.seed(123)
St = MC_gbm1(S0, mu, sigma, T, M, n)
# GBM Simulation 2
np.random.seed(123)
St2 = MC_gbm2(S0, mu, sigma, T, M, n)
print("Equal matrices?", np.allclose(St, St2)) # False
In this case, the difference comes from how you generate the t iteration in method 1: you use the previous value to compute the next one. Changing the corresponding line to
paths[t] = S0 * np.exp((mu- 0.5 * sigma ** 2) * dt + sigma * np.zeros_like(rand))
both functions will give the same result when the noise is turned off.
Finally, when you turn on the noise, the source of the difference is the extra step in method 2: the cumsum spoils the result. If you simply do
paths = S0 * np.exp((mu - 0.5 * sigma**2) * dt + sigma * dW)
in addition to the already discussed
paths[t] = S0 * np.exp((mu- 0.5 * sigma ** 2) * dt + sigma * rand)
you will get the same result:
# GBM Simulation 1
np.random.seed(123)
St = MC_gbm1(S0, mu, sigma, T, M, n)
# GBM Simulation 2
np.random.seed(123)
St2 = MC_gbm2(S0, mu, sigma, T, M, n)
print("Equal matrices?", np.allclose(St, St2)) # True