<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 データの取得方法
ここではQuandl.comからのデータを受け取っています。今回入手した日経平均株価は、
時間、開始値、最高値、最低値、終値のデータを入手していますが、古いデータは終値しかないようですので、終値を用います。
*** TODO いつからデータを入手することが最も効果的かを考える。（処理時間と制度に影響が出るため）
Step1:  抜けデータが目立ったため、週単位でのデータを入手します
Step2:  データの用い方
必要となるpythonパッケージのインポートを行っています。
*** TODO 実装はClojureで行いため、これに相当するパッケージを検索、作成を行う
Step3:  以下のグラフから、2000年ごろのデータからの推測でも十分に予測が行える可能性が伺えます。
Step4:  ARIMAモデルでモデル推定を行うための下準備として、株価の変化量を取得します。
Step5:  AICを求めてモデルの良さを計算しますが、やや時間（約三分）がかかってしまします。
（SARIMAモデルでこれを行うと、更に時間がかかります）
*** TODO 実行時間の計測と最適化・マシンスペックの向上と性能の関係の調査
Step6:  先程の実行結果から、AR=2, MA=2という値の場合が最も良いモデルになることがわかりました。
Step7:  比較のためSARIMAモデルではなく、ARIMAモデルでの推定を行ってみます。
こちらの実行はそれほど時間がかかりません。
Step8:  予測のブレがあまりないことが伺えます
Step9:  SARIMAモデルでの推定を行ってみます。
ARIMAモデルの実行がそれほど時間がかからなかったのに対して、SARIMAモデルはやや時間がかかること、Wariningが出ることが難点です。
Step10:  おおよそ見た限りではARIMAモデルと大差はないようですが、他の論文を読む限りではこちらの手法のほうが推測が上手く行くようです。
*** TODO データを少なくした場合の実行結果
Step11:  以下が予測を結合したもの
Step12:  青が実測値、赤が予測値です。それに近い値を計測できたのではないでしょうか？
後のために関数化しておく
<ASSISTANT_TASK:>
Python 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] # 最新のデータ１０件を表示
Explanation: データの取得方法
ここではQuandl.comからのデータを受け取っています。今回入手した日経平均株価は、
時間、開始値、最高値、最低値、終値のデータを入手していますが、古いデータは終値しかないようですので、終値を用います。
*** TODO いつからデータを入手することが最も効果的かを考える。（処理時間と制度に影響が出るため）
End of explanation
data_normal = data_normal.fillna(method='pad').resample('W-MON').fillna(method='pad')
data_normal[:5]
type(data_normal.index[0])
data_normal.index
Explanation: 抜けデータが目立ったため、週単位でのデータを入手します
End of explanation
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
Explanation: データの用い方
必要となるpythonパッケージのインポートを行っています。
*** TODO 実装はClojureで行いため、これに相当するパッケージを検索、作成を行う
End of explanation
plt.plot(data_normal)
Explanation: 以下のグラフから、2000年ごろのデータからの推測でも十分に予測が行える可能性が伺えます。
End of explanation
# 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)
Explanation: ARIMAモデルでモデル推定を行うための下準備として、株価の変化量を取得します。
End of explanation
# automatically ARIMA prediction function (using AIC)
resDiff = sm.tsa.arma_order_select_ic(diff, ic='aic', trend='nc')
# few Times ...(orz...)
Explanation: AICを求めてモデルの良さを計算しますが、やや時間（約三分）がかかってしまします。
（SARIMAモデルでこれを行うと、更に時間がかかります）
*** TODO 実行時間の計測と最適化・マシンスペックの向上と性能の関係の調査
End of explanation
resDiff 
# search min
resDiff['aic_min_order']
Explanation: 先程の実行結果から、AR=2, MA=2という値の場合が最も良いモデルになることがわかりました。
End of explanation
# 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
Explanation: 比較のためSARIMAモデルではなく、ARIMAモデルでの推定を行ってみます。
こちらの実行はそれほど時間がかかりません。
End of explanation
# 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
Explanation: 予測のブレがあまりないことが伺えます
End of explanation
# 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
Explanation: SARIMAモデルでの推定を行ってみます。
ARIMAモデルの実行がそれほど時間がかからなかったのに対して、SARIMAモデルはやや時間がかかること、Wariningが出ることが難点です。
End of explanation
# 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")
Explanation: おおよそ見た限りではARIMAモデルと大差はないようですが、他の論文を読む限りではこちらの手法のほうが推測が上手く行くようです。
*** TODO データを少なくした場合の実行結果
End of explanation
data_extra = pd.concat( [data_normal, pred[data_normal.index[-1] + 1:]] )
plt.plot(data_extra[:-150:-1])
Explanation: 以下が予測を結合したもの
End of explanation
# 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])
Explanation: 青が実測値、赤が予測値です。それに近い値を計測できたのではないでしょうか？
後のために関数化しておく
End of explanation
<END_TASK>