Impressum

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"))

mcsdaten =pd.read_csv("meininput/beispiel_stocks_5j.csv")

mcsdaten.head()

mcsdaten.tail()

#Es wird nur der Stock AAL ausgewählt
mcsdaten = mcsdaten[mcsdaten.Name == 'AAL']

mcsdaten.tail()

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

S0 = mcsdaten.close.iloc[-1]
S0

#Erst wird eine leere Matrix für die daily_returns erzeugt
price_list = np.zeros_like(daily_returns)
price_list[0] = S0
price_list

# 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()

close = mcsdaten.close
close = pd.DataFrame(close)
frames = [close, price_list]
monte_carlo_forecast = pd.concat(frames)

monte_carlo_forecast.head()

monte_carlo_forecast.tail()

monte_carlo = monte_carlo_forecast.iloc[:,:].values
import matplotlib.pyplot as plt
plt.figure(figsize=(17,8))
plt.plot(monte_carlo)
plt.show()

monte_carlo = monte_carlo_forecast.iloc[:,:].values
import matplotlib.pyplot as plt
plt.figure(figsize=(17,8))
plt.plot(monte_carlo)
plt.show()

#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)

mcsdaten = mcsdaten.dropna()

#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()

#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()

#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

#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())