Skip to content

Instantly share code, notes, and snippets.

@yildizumut yildizumut/GBM_sim.py Secret
Last active Aug 13, 2019

Embed
What would you like to do?
Full code
import pandas as pd
import numpy as np
import quandl
import matplotlib.pyplot as plt
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
stock_name = 'FSE/EON_X'
# Plot of E.ON(a big scale energy company in Europe)
#stock prices since beginning of 2019 (up to July)
prices = quandl.get(stock_name,
authtoken="your_authorization_token",
start_date = '2019-01-01', end_date = '2019-07-31'
).reset_index(drop = False)[['Date', 'Close']]
plt.figure(figsize = (15, 5))
plt.plot(prices['Date'], prices['Close'])
plt.xlabel('Days')
plt.ylabel('Stock Prices, €')
plt.show()
#--------------------------------------------------- GEOMETRIC BROWNIAN MOTION ------------------------------------------------
# Parameter Definitions
# So : initial stock price
# dt : time increment -> a day in our case
# T : length of the prediction time horizon(how many time points to predict, same unit with dt(days))
# N : number of time points in prediction the time horizon -> T/dt
# t : array for time points in the prediction time horizon [1, 2, 3, .. , N]
# mu : mean of historical daily returns
# sigma : standard deviation of historical daily returns
# b : array for brownian increments
# W : array for brownian path
start_date = '2019-07-01'
end_date = '2019-07-31'
pred_end_date = '2019-08-31'
# We get daily closing stock prices of E.ON for July 2019
S_eon = quandl.get(stock_name,
authtoken="your_authorization_token",
start_date = start_date, end_date = end_date
).reset_index(drop = False)[['Date', 'Close']]
print(S_eon.head())
print(S_eon.tail())
returns = (S_eon.loc[1:, 'Close'] - \
S_eon.shift(1).loc[1:, 'Close']) / \
S_eon.shift(1).loc[1:, 'Close']
print(returns.tolist())
# Parameter Assignments
So = S_eon.loc[S_eon.shape[0] - 1, "Close"]
dt = 1 # day # User input
n_of_wkdays = pd.date_range(start = pd.to_datetime(end_date,
format = "%Y-%m-%d") + pd.Timedelta('1 days'),
end = pd.to_datetime(pred_end_date,
format = "%Y-%m-%d")).to_series().map(lambda x:
1 if x.isoweekday() in range(1,6) else 0).sum()
T = n_of_wkdays # days # User input -> follows from pred_end_date
N = T / dt
t = np.arange(1, int(N) + 1)
mu = np.mean(returns)
sigma = np.std(returns)
scen_size = 50 # User input
b = {str(scen): np.random.normal(0, 1, int(N)) for scen in range(1, scen_size + 1)}
W = {str(scen): b[str(scen)].cumsum() for scen in range(1, scen_size + 1)}
# Calculating drift and diffusion components
drift = (mu - 0.5 * sigma**2) * t
print(drift)
diffusion = {str(scen): sigma * W[str(scen)] for scen in range(1, scen_size + 1)}
print(diffusion)
# Making the predictions
S = np.array([So * np.exp(drift + diffusion[str(scen)]) for scen in range(1, scen_size + 1)])
S = np.hstack((np.array([[So] for scen in range(scen_size)]), S)) # add So to the beginning series
print(S)
# Plotting the simulations
plt.figure(figsize = (20,10))
for i in range(scen_size):
plt.title("Daily Volatility: " + str(sigma))
plt.plot(pd.date_range(start = S_eon["Date"].max(),
end = pred_end_date, freq = 'D').map(lambda x:
x if x.isoweekday() in range(1, 6) else np.nan).dropna(), S[i, :])
plt.ylabel('Stock Prices, €')
plt.xlabel('Prediction Days')
plt.show()
# Dataframe format for predictions - first 10 scenarios only
Preds_df = pd.DataFrame(S.swapaxes(0, 1)[:, :10]).set_index(
pd.date_range(start = S_eon["Date"].max(),
end = pred_end_date, freq = 'D').map(lambda x:
x if x.isoweekday() in range(1, 6) else np.nan).dropna()
).reset_index(drop = False)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.