Solution

#Set a new parameter, runs which allows us to toggle the number of simulations
rf = .02
time_step = 1/252
price0 = 100
sigma = .08
runs = 100000
periods = 252

#Build a function to wrap everything in a function
def simulate_returns(price0, rf, time_step, periods, sigma):
    returns = pd.DataFrame(np.random.normal(0,1, periods * runs).reshape(periods, runs))
    returns = returns * sigma * timestep ** .5 + (rf - .5 * sigma ** 2) * timestep
    returns = np.exp(returns).cumprod()
    returns.index = returns.index + 1
    returns.loc[0] = 1
    returns = returns.sort_index()
    paths = returns * price0
    return paths
np.random.seed(1)
paths = simulate_returns(price0, rf, time_step, periods, sigma)
print(paths)

We will try to see how important using enough runs can be. We expected the stock price to go up linearly based on the risk-free rate. If we test with different numbers of test runs we see that as we add more and more runs we begin to converge to the correct linear shape for the average value of the stock across runs at each time step.

#Compare the accuracy with different number of simulations
mean_path = [paths.iloc[:,:x].mean(axis=1) for x in [100, 500, 1000, 10000, 100000]]
mean_path = pd.concat(mean_path, axis=1)
mean_path.columns = [100, 500, 1000, 10000, 100000]
mean_path.plot(kind='line')
plt.xlabel("Time Step")
plt.ylabel("Mean Stock Price")
plt.title("Monte Carlo Accuracy vs. # of Paths")
plt.show()

Another important part of this is that when we use this method we end up getting an option value similar to the black-scholes model. Recall what the price would be for a $110 strike price call.

#Recall the black-scholes formula
def black_scholes(S, X, sigma, rf, t):
    d1 = 1/(sigma*t**.5) * (np.log(S/X)  + (rf + sigma **2 /2) * t)
    d2 = d1 - sigma * t**.5
    return norm.cdf(d1) * S - norm.cdf(d2) * X * np.exp(-rf*t)
print(black_scholes(100, 110, .08, .02, 1))

Now, to value the call option, we need to take the value of a call option at the end of each run and then get the mean call option value at the end. Then, we can apply the discount factor to find the present value of the expected payoff on average based on all paths.

#Note the call option will be valued about the same
print(paths.iloc[-1].apply(lambda x: max(x-110, 0)).mean() / np.exp(.02 * 1))

The two match! You can also confirm that this is the case when using a strike price of $100.

#You can confirm with a different strike price of $100
print(black_scholes(100, 100, .08, .02, 1))
print(paths.iloc[-1].apply(lambda x: max(x-100, 0)).mean() / np.exp(.02))