<jupyter_start><jupyter_text>データの取得方法 ここではQuandl.comからのデータを受け取っています。今回入手した日経平均株価は、時間、開始値、最高値、最低値、終値のデータを入手していますが、古いデータは終値しかないようですので、終値を用います。*** TODO いつからデータを入手することが最も効果的かを考える。（処理時間と制度に影響が出るため）<jupyter_code>import quandl
data = quandl.get('NIKKEI/INDEX')
data[:5]
data_normal = (((data['Close Price']).to_frame())[-10000:-1])['Close Price']
data_normal[-10:-1] # 最新のデータ１０件を表示<jupyter_output><empty_output><jupyter_text>抜けデータが目立ったため、週単位でのデータを入手します<jupyter_code>data_normal = data_normal.fillna(method='pad').resample('W-MON').fillna(method='pad')
data_normal[:5]
type(data_normal.index[0])
data_normal.index<jupyter_output><empty_output><jupyter_text>データの用い方 必要となるpythonパッケージのインポートを行っています。*** TODO 実装はClojureで行いため、これに相当するパッケージを検索、作成を行う<jupyter_code>import numpy as np
import pandas as pd
from scipy import stats
from pandas.core import datetools
# grapgh plotting 
from matplotlib import pylab as plt
import seaborn as sns
%matplotlib inline

# settings graph size
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15,6

# model 
import statsmodels.api as sm<jupyter_output><empty_output><jupyter_text>以下のグラフから、2000年ごろのデータからの推測でも十分に予測が行える可能性が伺えます。<jupyter_code>plt.plot(data_normal)<jupyter_output><empty_output><jupyter_text>ARIMAモデルでモデル推定を行うための下準備として、株価の変化量を取得します。<jupyter_code># ARIMA model prediction ... (This is self thought (not automatically))
diff = data_normal - data_normal.shift()
diff = diff.dropna()
diff.head()
# difference plot 
plt.plot(diff)<jupyter_output><empty_output><jupyter_text>AICを求めてモデルの良さを計算しますが、やや時間（約三分）がかかってしまします。（SARIMAモデルでこれを行うと、更に時間がかかります）*** TODO 実行時間の計測と最適化・マシンスペックの向上と性能の関係の調査<jupyter_code># automatically ARIMA prediction function (using AIC)
resDiff = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')
# few Times ...(orz...)<jupyter_output><empty_output><jupyter_text>先程の実行結果から、AR=2, MA=2という値の場合が最も良いモデルになることがわかりました。<jupyter_code>resDiff 
# search min
resDiff['aic_min_order']<jupyter_output><empty_output><jupyter_text>比較のためSARIMAモデルではなく、ARIMAモデルでの推定を行ってみます。 こちらの実行はそれほど時間がかかりません。<jupyter_code># we found x = x, y= y autopmatically
from statsmodels.tsa.arima_model import ARIMA
ARIMAx_1_y = ARIMA(data_normal, 
                  order=(resDiff['aic_min_order'][0], 1,
                  resDiff['aic_min_order'][1])).fit(dist=False)
# AR = resDiff[...][0] / I = 1 / MA = resDiff[...][1]
ARIMAx_1_y.params<jupyter_output><empty_output><jupyter_text>予測のブレがあまりないことが伺えます<jupyter_code># check Residual error (... I think this is "White noise")
# this is not Arima ... (Periodicity remained)
resid = ARIMAx_1_y.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)

# ok? 
# We test SARIMA_model<jupyter_output><empty_output><jupyter_text>SARIMAモデルでの推定を行ってみます。 ARIMAモデルの実行がそれほど時間がかからなかったのに対して、SARIMAモデルはやや時間がかかること、Wariningが出ることが難点です。<jupyter_code># predict SARIMA model  by  myself (not automatically) 
import statsmodels.api as sm

SARIMAx_1_y_111 = sm.tsa.SARIMAX(data_normal, 
                                 order=(2,1,2),seasonal_order=(1,1,1,12))
SARIMAx_1_y_111 = SARIMAx_1_y_111.fit()
# order ... from ARIMA model // seasonal_order ... 1 1 1 ... ?  
print(SARIMAx_1_y_111.summary())

# maybe use "Box-Jenkins method" ...
# https://github.com/statsmodels/statsmodels/issues/3620 for error<jupyter_output>Statespace Model Results                                 
==========================================================================================
Dep. Variable:                        Close Price   No. Observations:                 2114
Model:             SARIMAX(2, 1, 2)x(1, 1, 1, 12)   Log Likelihood              -15992.035
Date:                            Sun, 06 Aug 2017   AIC                          31998.069
Time:                                    19:01:40   BIC                          32037.663
Sample:                                02-07-1977   HQIC                         32012.567
                                     - 08-07-2017                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
-----------------------[...]<jupyter_text>おおよそ見た限りではARIMAモデルと大差はないようですが、他の論文を読む限りではこちらの手法のほうが推測が上手く行くようです。*** TODO データを少なくした場合の実行結果<jupyter_code># check Residual error
residSARIMA = SARIMAx_1_y_111.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residSARIMA.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residSARIMA, lags=40, ax=ax2)
# prediction
pred = SARIMAx_1_y_111.predict(start = 1, end = '2018-01-15')
# (print(SARIMAx_1_y_111.__doc__))

# 本来は未来（インデクスの外）まで予測ができるはずなのですが、
# 何故かエラーが出てしまうので、既存のデータ部分だけ予測します

# TODO エラーの原因特定
# plot real data and predict data 
plt.plot(data_normal[:-150:-1])
plt.plot(pred[:-150:-1], "r")<jupyter_output><empty_output><jupyter_text>以下が予測を結合したもの<jupyter_code>data_extra = pd.concat( [data_normal, pred[data_normal.index[-1] + 1:]] )
plt.plot(data_extra[:-150:-1])<jupyter_output><empty_output><jupyter_text>青が実測値、赤が予測値です。それに近い値を計測できたのではないでしょうか？  後のために関数化しておく<jupyter_code># require 
import quandl
import numpy as np
import pandad as pd 
from scipy import stats
from pandas.coda import datatools
import statsmodels.api as sm
def get_data(quandl_name):
    data = quandl.get(quandl_name)
    return data

def set_data(data):
    data_normal = (((data['Close Price']).to_frame())[-10000:-1])['Close Price']
    data_normal = data_normal.fillna(method='pad').resample('W-MON').fillna(method='pad')
    return data_normal

def sarima(quandl_name):
    data_normal = set_data(get_data(quandl_name))
    diff = (data_normal - (data_normal.shift())).dropna()
    resDiff = aic(diff)['aic_min_order']
    ar = resDiff[0]
    ma = resDiff[1]
    SARIMAx_1_y_111 =  \
    sm.tsa.SARIMAX(data_normal, order=(int(ar),1, int(ma)),seasonal_order=(1,1,1,12))
    return  SARIMAx_1_y_111

def pred_data(SARIMAx_1_y_111, predict_date):
    SARIMAx_1_y_111 = SARIMAx_1_y_111.fit()
    print(SARIMAx_1_y_111.summary())
    pred = SARIMAx_1_y_111.predict(start = 1, end = predict_date)
    return pd.concat( [data_normal, pred[data_normal.index[-1] + 1:]] )

def aic (diff):
    return (sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')) 


# 以上をまとめたもの
def predict_data(quandl_name, predict_data):
    sarima_model = sarima(quandl_name)
    return pred_data(sarima_model, predict_data)
predict_res = predict_data('NIKKEI/INDEX','2018-01-15')
plt.plot(predict_res[:-150:-1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1488673139158576, 'fraction_numerical': 0.03957046190055899, 'mean_word_length': 4.2179585571757485, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 53, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '467589', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3239, 'n_tokens_neox': 2969, 'n_words': 379, 'path': 'value-tracker/report/python/.ipynb_checkpoints/report-checkpoint.ipynb', 'repo_name': 'elect000/Journal', 'size': '6798'}