In [1]:
import numpy as np # Import des Basis-Pakets
import pandas as pd # Import zur Basis-Datenverarbeitung, CSV-Datei lesen und schreiben etc.(z. B. per pd.read_csv)
# Input-Daten sind in dem Ordner "meininput/" verfügbar.
# Der nachfolgende Quellcode listet alle in dem Ordner enthaltenen Dateien.
import os
print(os.listdir("meininput"))
In [3]:
mcsdaten =pd.read_csv("meininput/beispiel_stocks_5j.csv")
In [4]:
mcsdaten.head()
Out[4]:
In [5]:
mcsdaten.tail()
Out[5]:
In [6]:
#Es wird nur der Stock AAL ausgewählt
mcsdaten = mcsdaten[mcsdaten.Name == 'AAL']
In [7]:
mcsdaten.tail()
Out[7]:
In [8]:
from scipy.stats import norm
log_returns = np.log(1 + mcsdaten.close.pct_change())
u = log_returns.mean() #Mean des log-Returns
var = log_returns.var() #Varianz des log-Returns
drift = u - (0.5 * var) #drift bzw. Trend des log-Returns
stdev = log_returns.std() #Standardabweichung des log-Returns
t_intervals = 250 #Forecast von 250 Zeitpunkten
iterations = 10 #10 unterschiedliche Forecasts (als loop)
daily_returns = np.exp(drift + stdev * norm.ppf(np.random.rand(t_intervals, iterations)))
#daily_returns als "noise", bei Multiplikation mit t time price, erhält man t+1 time price
In [9]:
S0 = mcsdaten.close.iloc[-1]
S0
Out[9]:
In [10]:
#Erst wird eine leere Matrix für die daily_returns erzeugt
price_list = np.zeros_like(daily_returns)
price_list[0] = S0
price_list
Out[10]:
In [11]:
# Mit einem einfachen Loop wird ein Forecast für die nächsten 250 Tage performt.
for t in range(1, t_intervals):
price_list[t] = price_list[t - 1] * daily_returns[t]
price_list = pd.DataFrame(price_list)
price_list['close'] = price_list[0]
price_list.head()
Out[11]:
In [12]:
close = mcsdaten.close
close = pd.DataFrame(close)
frames = [close, price_list]
monte_carlo_forecast = pd.concat(frames)
In [13]:
monte_carlo_forecast.head()
Out[13]:
In [14]:
monte_carlo_forecast.tail()
Out[14]:
In [15]:
monte_carlo = monte_carlo_forecast.iloc[:,:].values
import matplotlib.pyplot as plt
plt.figure(figsize=(17,8))
plt.plot(monte_carlo)
plt.show()
In [16]:
monte_carlo = monte_carlo_forecast.iloc[:,:].values
import matplotlib.pyplot as plt
plt.figure(figsize=(17,8))
plt.plot(monte_carlo)
plt.show()
In [19]:
#Nachfolgend die Verteilung der log-Returns
import plotly.plotly as py
from plotly.offline import init_notebook_mode, iplot
init_notebook_mode(connected=True)
import plotly.graph_objs as go
trace = go.Histogram(x=log_returns,opacity=0.85,name = "Logarithmic Return", marker=dict(color='rgba(0, 0, 255, 0.8)'))
info = [trace]
layout = go.Layout(barmode='overlay',
title='Verteilung der Logarithmic Returns',
xaxis=dict(title='Logarithmic Return'),
yaxis=dict( title='Verteilung'),
)
fig = go.Figure(data=info, layout=layout)
iplot(fig)
In [26]:
mcsdaten = mcsdaten.dropna()
In [28]:
#Nun wird die Seasonalität untersucht
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(mcsdaten.log_return, freq = 260) #mit 260 Arbeitstagen pro Jahr
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.figure(figsize=(12,6))
plt.subplot(411)
plt.plot(log_returns, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonalität')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuale')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
In [29]:
#Nun wird die serielle Korrelation untersucht
import statsmodels.api as sm
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.stattools import pacf
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(mcsdaten.log_return, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(mcsdaten.log_return, lags=40, ax=ax2)
plt.show()
In [30]:
#Nun wird das Auto Regressive Moving Average Modell zugrunde gelegt
#...mit den Parametern Ar(p) and Ma(q)
from statsmodels.tsa.stattools import ARMA
def best_AR_MA_checker(df,lower,upper):
from statsmodels.tsa.stattools import ARMA
from statsmodels.tsa.stattools import adfuller
arg=np.arange(lower,upper)
arg1=np.arange(lower,upper)
best_param_i=0
best_param_j=0
temp=12000000
rs=99
for i in arg:
for j in arg1:
model=ARMA(df, order=(i,0,j))
result=model.fit(disp=0)
resid=adfuller(result.resid)
if (result.aic<temp and adfuller(result.resid)[1]<0.05):
temp=result.aic
best_param_i=i
best_param_j=j
rs=resid[1]
print ("AR: %d, MA: %d, AIC: %d; resid stationarity check: %d"%(i,j,result.aic,resid[1]))
print("the following function prints AIC criteria and finds the paramters for minimum AIC criteria")
print("best AR: %d, best MA: %d, best AIC: %d; resid stationarity check:%d"%(best_param_i, best_param_j, temp, rs))
best_AR_MA_checker(mcsdaten.log_return,0,3) #For each parameter I want to try from 0 to 2
In [31]:
#Nun wird arma(1,0) untersucht
from statsmodels.tsa.stattools import ARMA
model=ARMA(mcsdaten.log_return, order=(1,0))
res=model.fit(disp=0)
print (res.summary())
In [ ]: