時系列予測をGBDTで実装してみた

* 製造業出身のデータサイエンティストがお送りする記事
* 今回は時系列データを回帰モデル（GBDT）を用いて予測してみました。

##はじめに
過去に時系列解析や回帰モデルの手法を整理しておりますので、興味ある方はそちらも参照して頂けますと幸いです。

* [時系列解析](https://qiita.com/DS27/items/1e998a58488e76bfcbdc)
* [状態空間モデル](https://qiita.com/DS27/items/cd0e11c2b9eb097c606b)
* [回帰モデル](https://qiita.com/DS27/items/aa3f6d0f03a8053e5810)

##GBDTの時系列予測

pythonのコードは下記の通りです。

```php
# 必要なライブラリーのインポート
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from matplotlib import pylab as plt
%matplotlib inline

# 統計モデル
import statsmodels.api as sm

# GBDT
from sklearn.ensemble import GradientBoostingRegressor

# グラフを横長にする
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

# https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html
df = pd.read_csv('AirPassengers.csv')

# float型に変換
df['#Passengers'] = df['#Passengers'].astype('float64')
df = df.rename(columns={'#Passengers': 'Passengers'})

# datetime型にしてインデックスにする
df.Month = pd.to_datetime(df.Month)
df = df.set_index("Month")

# データの中身を確認
df.head()
```
![スクリーンショット 2020-12-05 12.41.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/e5bf2b87-d8c3-49ef-66a9-7e697e6c9b7f.png)

次にコレログラムを作成します。

```php
#  自己相関のグラフ
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(df["Passengers"], lags=30)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/44814e46-dae8-c2da-9b2a-5090b4539e05.png)

```php
# 偏自己相関を可視化する
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_pacf(df["Passengers"], lags=20)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/d704d069-e244-7fdd-a1e2-02bc06b3f133.png)

今回のデータでは、偏自己相関のグラフを見ると、12ヶ月周期で相関があることが分かるかと思います。つまり、季節的な周期変動があることが分かります。

次に過去12ヶ月の履歴を作成します。

```php
for i in range(1, 13):
    df['shift%s'%i] = df['Passengers'].shift(i)

pd.concat([df.head(13), df.tail(3)], axis=0, sort=False)
```

![スクリーンショット 2021-01-11 13.30.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/419cb658-1e83-de4f-ac88-8f6d54c226d4.png)

次に時系列データでよく使用する差分列を作成します。

```php
df['deriv1'] = df['shift1'].diff(1)
df[['Passengers', 'deriv1']].head()
```

![スクリーンショット 2021-01-11 13.31.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/33438448-81da-4dfe-c808-49bfddc33f56.png)

次は2回差分列を作成します。

```php
df['deriv2'] = df['shift1'].diff(1).diff(1)
df[['Passengers', 'deriv2']].head()
```

![スクリーンショット 2021-01-11 13.32.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/b57bd643-de18-dfd3-bad3-49409dcf9815.png)

最後に統計量も説明変数に追加します。

```php
df['mean'] = df['shift1'].rolling(12).mean()
df['median'] = df['shift1'].rolling(12).median()
df['max'] = df['shift1'].rolling(12).max()
df['min'] = df['shift1'].rolling(12).min()
df[['Passengers', 'mean', 'median', 'max', 'min']][12:24]
```

![スクリーンショット 2021-01-11 13.33.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/dbea170c-2b2e-bed4-11de-f461c1f7701d.png)

これからGBDTで予測をしていきます。

```php
# 欠損値データの削除
df = df.dropna()
df.head()

x = df.drop('Passengers', axis=1)
y = df['Passengers']

# 学習データと評価データを作成
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=1)

#データを標準化
sc = StandardScaler()
sc.fit(x_train) #学習用データで標準化
x_train_std = sc.transform(x_train)
x_test_std = sc.transform(x_test)

# モデルの学習
GBDT = GradientBoostingRegressor()
GBDT.fit(x_train_std, y_train)

# 予測　
y_pred = GBDT.predict(x_test_std)

y_ = np.concatenate([np.array([None for i in range(len(y_train))]), y_pred])
y_ = pd.DataFrame(y_, index=df.index)

plt.figure(figsize=(10,5))
plt.plot(y, label='original')
plt.plot(y_, '--', label='predict')
plt.legend()
```

![結果.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/743453/17f624eb-e9ba-b1cc-d001-1209df278b52.png)

## さいごに
最後まで読んで頂き、ありがとうございました。
今回は時系列データを回帰モデルを用いて予測してみました。
回帰モデルを使用する際は、特徴量作成と選択が重要ですね。

訂正要望がありましたら、ご連絡頂けますと幸いです。