時系列解析ライブラリーの比較 statsmodel,fbprophet,tensorflow　その2

# はじめに
　前回の投稿の続きとして、tensorflowを用いて時系列予測に取り組む。対象は、前回と同様に東京電力の公開する電力需要データとする。

# モデル化
　前回のstatsmodelを用いたモデル化の際に得られた知見をもとにモデル化する。以下に、電力需要データの特徴とともにモデル化の際の考慮事項を記載する。

1.短期周期性
　過去2期との自己相関がうかがえることからこれらをモデル化に組み込むことにする。イメージとしては、

　　　電力需要t ＝ 電力需要t-1 × Weight + 電力需要t-2 × Weight ⁺ Bias

ということで、変数ふたつのニューラルネットワークとする。隠れ層はあまり複雑にせず1階とする。狙いとしては、過去2期のデータから本質的な特徴を抽出し、それを予測に使うというものである。隠れ層のノード数は、週平均にまとめてあるとはいえ、変化がそれなりに激しいことから複雑度に合わせて100とする。活性化関数は、正負の値をうまく表現するtanhを使うことで2期の値が上げ傾向なのか下げ傾向なのかを表現する。

2.季節周期性
　以下のグラフは、電力需要を月別平均にしてプロットしたものである。

![fig52.png](https://qiita-image-store.s3.amazonaws.com/0/163753/486dce90-8a66-9154-489b-e8718cd710a5.png)

　これから、1年周期の存在が明らかである。つまり、1年のうちどの月に属するのかがわかれば、それによって平均的な電力需要量が推測できる。この月をカテゴリーデータとしてモデル化に取り込むことでこの月の影響を電力需要に反映する。隠れ層であるが、短期自己相関に合わせて1層、ノード数100としてネットワーク化を図る。

　以上をコードにしたものが以下である。

```
x1 = tf.placeholder(tf.float32,[None,2])
x2 = tf.placeholder(tf.float32,[None,12])

w11 = tf.Variable(tf.truncated_normal([2,100]))
b11 = tf.Variable(tf.zeros([100]))
w12 = tf.Variable(tf.truncated_normal([12,100]))
b12 = tf.Variable(tf.zeros([100]))

w21 = tf.Variable(tf.truncated_normal([100,1]))
b21 = tf.Variable(tf.zeros([1]))
w22 = tf.Variable(tf.truncated_normal([100,1]))
b22 = tf.Variable(tf.zeros([1]))

y01 = tf.nn.tanh(tf.matmul(x1 ,w11) + b11)
y02 = tf.nn.tanh(tf.matmul(x2 ,w12) + b12)

y11 = tf.matmul(y01,w21) + b21
y12 = tf.matmul(y02,w22) + b22

y1 = y11 + y12

y = tf.placeholder(tf.float32,[None,1])

loss = tf.reduce_sum(tf.square(y1-y))

train_step = tf.train.AdamOptimizer().minimize(loss) 
```
 なお、短期自己相関と季節周期性は、独立して予測結果に影響していると仮定して、短期自己相関による結果（y11）と季節周期性による結果（y12）を別にしている。

# データの生成

　先のモデルに合わせたデータの生成のコードを以下に掲載する。

```
url = "http://www.tepco.co.jp/forecast/html/images/"
file_list =["juyo-2016.csv","juyo-2017.csv","juyo-2018.csv"]

dat = pd.DataFrame([[np.nan,np.nan,np.nan]])

for i in file_list :  #東京電力から電力需要値を取得する
    res = requests.get(url+i)
    docs = res.text.split('\n')[3:]
    tmp = pd.DataFrame([x.split(',') for x in docs])
    tmp.dropna(inplace=True)
    tmp[2] = [int(x.replace("\r","")) for x in tmp[2]]
    dat = pd.concat([dat,tmp])

dat.dropna(inplace=True)

dat.reset_index(drop=True,inplace=True)
dat.columns = ['Date','Time','実績']

dat['Date_Time'] = dat['Date'] + ' ' + dat['Time']  #日付をindexに設定
dat.index = pd.to_datetime(dat['Date_Time'])

dat.drop('Date_Time',axis=1,inplace=True)

dat_w = dat.resample('W').mean()  #週平均のデータに変換

dat_w['月'] = [x.month for x in dat_w.index]
time_df = dat_w.copy()
tmp = pd.get_dummies(time_df['月'])
tmp.columns = ['Jan','Feb','Mar','Apl','May','Jun','Jly','Aug','Sep','Oct','Nov','Dec']
time_df = pd.concat([time_df,tmp],axis=1)
time_df['t-1'] = time_df.shift(1)['実績']
time_df['t-2'] = time_df.shift(2)['実績']

time_df.dropna(inplace=True)

X1 = time_df[['t-2','t-1']].as_matrix()
X2 = time_df[['Jan','Feb','Mar','Apl','May','Jun','Jly','Aug','Sep','Oct','Nov','Dec']].as_matrix()
Y = time_df['実績'].as_matrix().reshape(len(time_df),1)

X1_train = X1[:70]
X2_train = X2[:70]
Y_train = Y[:70]
```
tensorflowに用いるデータはarrayなのだが、arrayの状態のままで加工するのは不便であるため、DataFrameを活用してデータを生成している。時系列データは、並びに意味があるため、シャッフルせずに70件目までを学習用データとして用いることにした。

# 学習の実施

```
sess = tf.Session()

sess.run(tf.global_variables_initializer())

print('Start',datetime.datetime.now())

i=0
loss_val = []
length = len(X1_train)
for j in range(12000):
    
    idx = 0
    for k in range(10) :
        if length - idx >=26 :
            X1_train_batch = X1_train[idx:idx+26]
            X2_train_batch = X2_train[idx:idx+26]
            Y_train_batch  = Y_train[idx:idx+26]
            idx = idx+26
        else :
            X1_train_batch = X1_train[idx:]
            X2_train_batch = X2_train[idx:]
            Y_train_batch  = Y_train[idx:]
            idx = 0
            
        sess.run(train_step , 
                 feed_dict={x1:X1_train_batch ,x2:X2_train_batch, y:Y_train_batch}) 

    i=i+1
    if i % 10 == 0 :
        tmp = sess.run(loss,
                       feed_dict={x1:X1_train_batch ,x2:X2_train_batch, y:Y_train_batch})    　
        loss_val.append([i,tmp])
print('End',datetime.datetime.now())
```
学習はミニバッチを応用している。データの周期性からすると半年（26週）で上げ下げを繰り返しているようなので、これを基準にデータを切り分けて学習させる。これによって季節周期性がよりうまく学習できると考えた。

# 結果
　まず、損失関数の収束状況を示す。
![fig53.png](https://qiita-image-store.s3.amazonaws.com/0/163753/f86a8756-1b74-6a57-10c8-7969f3b0cd3b.png)

　順調に減少しており、学習はうまくいったようである。
　続いてモデルの予測結果を見たい。対象は全期間とする。これによって学習に用いた期間とそれ以後の期間の予測精度を確認する。

![fig51.png](https://qiita-image-store.s3.amazonaws.com/0/163753/2529f3e4-5a0a-9025-5e63-c035588104ab.png)
図中の赤線は学習に用いた期間と試験に用いた期間の境界を示す。これからすると、学習と試験でそれほどの予測精度の相違はなく、汎化性の高いモデルになっていると言えそうである。

# まとめ
　tensorflowの場合、モデルの作り方はかなりの自由度がある。それゆえに利用者の思いついたアイデアを任意にモデルとして表現することが可能である。一方、自由度が高いがゆえに精度を究極的にまで突き詰めようとした場合、考慮すべき事項が多数にわたるため、苦労しそうである。
　今回は2期までを特徴量としたが、さらに広げることもできる。また、移動平均はあえて使わないモデルにしているが、織り込むことも可能である。また、一般的なニューラルネットと同様にミニバッチのサイズやエポック数などの配慮も性能に影響することが考えられる。
　さらにいうと、RNNやLSTMなどの時系列モデルに適用したレイヤーを組み込むことで精度向上も可能であろう。ただ、RNNの組み込みは思いのほか簡単ではなさそうだ。当初、先のコードにRNNセルを表す関数を追記すれば実装できる見込んでいたが、うまく機能せず、改めてRNN実装に関する知識を積んでから取り組んだほうがいいことを認識させられた。近いうちに習作ができたら報告したい。