回帰モデルの比較 - ARMA vs. Random Forest Regression


（追記）
本記事後半部，Random Rorest Regression に関する訂正記事を投稿しました．こちらも参照ください．[時系列データ分析の処理でやってはいけないこと（反省を含めて）](http://qiita.com/TomokIshii/items/ac7bde63f2c0e0de47b3)

Pythonで時系列データ分析を行う場合，統計モデリングのStatsModelsが備えるTimes series analysis (tsa) モジュールを使うのが第1の選択肢になると考えられる．また，機械学習ライブラリ  Scikit-learn には別のコンセプトで回帰を行う，決定木回帰(Decision Tree Regression)やRandom Forest Regressionがあるので，これも試してみたいと考え，比較検討した．

（プログラミング環境：Python 2.7.11, Pandas 0.18.0, StatsModels 0.6.1, Scikit-learn 0.17.1 になります．）

## 伝統的な自己回帰モデルのグループ

時系列データ分析の文献を開くと多くの場合，次の順番で説明が出てくる．

| 略称 | 説明  |
|:-----:|:-----|
| **AR** | 自己回帰(Auto Regression)モデル |
| **MA** | 移動平均(Moving Average)モデル |
| **ARMA** | 自己回帰移動平均モデル  |
| **ARIMA** | 自己回帰和分移動平均モデル |

各モデルの詳細はWikipedia等を参照していただくとして，基本的に**ARMA**は，**AR** , **MA**を含有するモデルであり， **ARIMA** は，**ARMA** を含むモデルなので，ライブラリとしては **ARIMA** がサポートされていれば，上記４つのモデルはすべて対応可能となる．但し，StatsModelsではAPIとして**AR**, **ARMA**, **ARIMA** が用意されている．

http://statsmodels.sourceforge.net/stable/tsa.html

今回，分析を行うデータとして，"Nile" を取り上げた．

> Nile River flows at Ashwan 1871-1970
> This dataset contains measurements on the annual flow of the Nile as measured at Ashwan for 100 years from 1871-1970. There is an apparent changepoint near 1898.

アフリカ，ナイル川流量の各年ごとのデータで，StatsModelライブラリ内に用意されているサンプルデータである．

![Nile_1.png](https://qiita-image-store.s3.amazonaws.com/0/74152/a012f36f-bf1f-470c-6f29-526bacf8f9d9.png)

自己回帰モデルでは多くの場合「定常性」がモデル適用の前提条件になっているが，上のプロットでは単調増加，単調減少の様子が見られないので，**ARMA**モデルを当てはめることにする．


### ARMAモデルの次数選択

およそ100年間の時系列データとなっているので，フィット用（訓練用）のデータとして前半70年分，モデルテスト用のデータとして後半30年分を割り当てることとする．ARMAモデルの次数(p, q)を選ぶ必要があるが，まずコレログラムとも呼ばれるACFプロット(Autocorrelation function plot), PACFプロット(Partial autocorrelation plot)を作図してみる．

```py
def load_data():
    df_read = sm.datasets.nile.load_pandas().data
    s_date = pd.Series(
        [pd.to_datetime(str(int(y_str))) for y_str in df_read['year']]
    )
    df = df_read.set_index(s_date)
    df = df.drop('year', axis=1)
    
    return df

df_nile = load_data()

# Plot Time Series Data
ax = df_nile['volume'].plot(figsize=(8,4), grid=True)
ax.set_xlabel('year')
ax.set_ylabel('volume x 10^8 m^3')
plt.show()

# Data Split (70: 30)
# ACF, PACF
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_train.values, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_train.values, lags=40, ax=ax2)
```

**Fig. ACF and PACF plot of "Nile" data**
![Nile_2.png](https://qiita-image-store.s3.amazonaws.com/0/74152/98f04c1b-ec1d-f685-d810-53b4c4b7e709.png)

ARMA(p, q) で，ACFプロット（上）より q=[0, 1, 2, 3] が選択肢に入り，PACFプロット（下）より p = [0, 1] が選択肢に入ることが（なんとなく）分かる．一応，ARMA(p=1, q=3) を視野に入れつつ次数選択のユーティリティで情報を出力させると以下のようになった．

```py
info_criteria = sm.tsa.stattools.arma_order_select_ic(
                    df_train.values, ic=['aic', 'bic']
                )
print(info_criteria.aic_min_order)
print(info_criteria.bic_min_order)

ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
>>> (1, 1)
>>> (1, 0)
```

収束性に関するWarningが出ていることに注意．Warningをふまえた上で，AIC（モデル選択基準値）からはARMA(1,  1), BIC（モデル選択基準値）からはAMMA(1, 0) が良さそうとの情報を得た．よってこの２つのモデルをデータにフィットさせてみる．

```py
arma_10 = sm.tsa.ARMA(df_train, (1, 0)).fit()
arma_11 = sm.tsa.ARMA(df_train, (1, 1)).fit()
```

これは，特にError, Warningもなく処理された．戻り値 arma_10, arma_11 はフィット後のモデル・オブジェクト(ARMAResultsクラス）である．一応，より高次のモデルARMA(1, 2), ARMA(1, 3) も試してみた．

```py
arma_12 = sm.tsa.ARMA(df_train, (1, 2)).fit()
arma_13 = sm.tsa.ARMA(df_train, (1, 3)).fit()

（中略，Traceback 情報など．）

ValueError: The computed initial AR coefficients are not stationary
You should induce stationarity, choose a different model order, or you can
pass your own start_params.

```

上のようなValueErrorが発生した．（エラーメッセージは，定常性(stationary)について言及していますが...）パラメータ算定時における収束性に問題がありARMA(1, 2)，ARMA(1, 3)への当てはめは難しい状況である．この先，この原因を深く追及してもいいが，今回は，ARMA vs. Random Forest Regressor がテーマであるので，フィットができたARMA(1, 0)とARMA(1, 1) を用いて検討を進める．

## ARMAモデルの検証
まず，AICの値を確認する．

```py
print('ARMA(1,0): AIC = %.2f' % arma_10.aic)
print('ARMA(1,1): AIC = %.2f' % arma_11.aic)
>>> ARMA(1,0): AIC = 910.94
>>> ARMA(1,1): AIC = 908.92
```

AIC値は，ARMA(1, 1)モデルの方が小さい．

次に訓練データ区間（70年分）のモデル推定値とその後のテスト区間（30年）のモデル推定値を算出する．

```py
# in-sample predict
arma_10_inpred = arma_10.predict(start='1871-01-01', end='1940-01-01')
arma_11_inpred = arma_11.predict(start='1871-01-01', end='1940-01-01')

# out-of-sample predict
arma_10_outpred = arma_10.predict(start='1941-01-01', end='1970-01-01')
arma_11_outpred = arma_11.predict(start='1941-01-01', end='1970-01-01')
```

AICの算定，推定値の算出ともにARMAモデル結果クラス(ARMAResultsクラス)のメソッドで求めることができるような実装となっている．元データの値，推定値をまとめてプロットする．

```py
# plot data and predicted values
def plot_ARMA_results(origdata, pred10in, pred10out, pred11in, pred11out):
    px = origdata.index
    py1 = origdata.values
    plt.plot(px, py1, 'b:', label='orig data')
    
    px_in = pred10in.index
    plt.plot(px_in, pred10in.values, 'g')
    plt.plot(px_in, pred11in.values, 'c')
    
    px_out = pred10out.index
    plt.plot(px_out, pred10out.values, 'g', label='ARMA(1,0)')
    plt.plot(px_out, pred11out.values, 'c', label='ARMA(1,1)')
    
    plt.legend()
    plt.grid(True)
    plt.show()

plot_ARMA_results(df_nile, arma_10_inpred, arma_10_outpred, 
    arma_11_inpred, arma_11_outpred)

```
**Fig. Nile data - original data and ARMA estimated data**
![Nile_3.png](https://qiita-image-store.s3.amazonaws.com/0/74152/8388894d-4633-8f50-9ebe-e91d73ae427d.png)

1870年から1940年までが訓練区間 (in-sample)，それ以降がテスト区間(out-of-sample)である．少し分かりにくいが，なんとなくARMA(1,1)の方が，実データに近いように見える．後でモデル比較のため，残差からMSE(Mean Squared Error)を求めておく．

```py
# Residue (mse for train) 
arma_10_mse_tr = np.array([r ** 2 for r in arma_10.resid]).mean()
arma_11_mse_tr = np.array([r ** 2 for r in arma_11.resid]).mean()

# Residue (mse for test)
arma_10_mse_te = np.array([(df_test.values[i] - arma_10_outpred[i]) **2 
                            for i in range(30)]).mean()
arma_11_mse_te = np.array([(df_test.values[i] - arma_11_outpred[i]) **2 
                            for i in range(30)]).mean()

```

## Random Forest 回帰

Random Forest回帰のベースとなっているのは，決定木回帰(Decision Tree Regression)である．決定木のアルゴリズム自体，大雑把でも理解しておく必要があると思われるが，本記事では説明を省略し，Scikit-learnの機能をブラックボックスとして扱うことにする．Scikit-learnのドキュメントにはExampleがあるが，それを試した結果が下図である．

**Fig. Decision Tree Regression Example**
![tree_regression_s.png](https://qiita-image-store.s3.amazonaws.com/0/74152/e9f130ef-06e3-a002-b92e-0833e46b0ded.png)

階段状の直線でフィットすることが特徴となっている．今回は単変量の時系列データであるが，いくつかの過去のデータ使って現在の値を推定するモデルとしてみた．具体的には，
` Volume_current ~ Volume_Lag1 + Volume_Lag2 + Volume_Lag3`
というモデルである．訓練データとテストデータを前回と同じ(70, 30)にしたかったが，Lag値を計算する過程でNaNがでるので，初めの部分をdropna()して，(67, 30) のデータ長とした．

まず，データの前処理を行う．

```py
df_nile['lag1'] = df_nile['volume'].shift(1) 
df_nile['lag2'] = df_nile['volume'].shift(2)
df_nile['lag3'] = df_nile['volume'].shift(3)

df_nile = df_nile.dropna()

X_train = df_nile[['lag1', 'lag2', 'lag3']][:67].values
X_test = df_nile[['lag1', 'lag2', 'lag3']][67:].values

y_train = df_nile['volume'][:67].values
y_test = df_nile['volume'][67:].values
```

Scikit-learnのRandomForestRegressorを使用．

```py
from sklearn.ensemble import RandomForestRegressor
r_forest = RandomForestRegressor(
            n_estimators=100,
            criterion='mse',
            random_state=1,
            n_jobs=-1
)
r_forest.fit(X_train, y_train)
y_train_pred = r_forest.predict(X_train)
y_test_pred = r_forest.predict(X_test)
```

推定値は，下図のようになった．

**Fig. Nile data - original data and Random Forest Regression results**
![Nile_4.png](https://qiita-image-store.s3.amazonaws.com/0/74152/c0f1e48a-a2c0-99fc-0a03-e3dd1eec4479.png)

最後にMSEを算定する．

```py
# check residue (mse)
train_resid = y_train - y_train_pred
RFR_mse_train = np.array([r ** 2 for r in train_resid]).mean()
test_resid = y_test - y_test_pred
RFR_mse_test = np.array([r ** 2 for r in test_resid]).mean()
```

## モデル精度比較
各モデルのMSE(Mean Squared Error)を比較する．

|    Model   |  MSE of in-samle (train)  | MSE of out-of-sample (test) |
|:----------:|:-----:|:-----:|
| ARMA(1,0) | 24111.0 | 18684.4 |
| ARMA(1,1) | 22757.3  | 16625.8 | 
| Random F. Regressor | 3470.1 | 15400.3 |

(MSE : smaller ... better)

次数の異なるARMAモデル同士の比較では，ARMA(1,1) の方が訓練区間，テスト区間ともMSEが小さく精度が良いことが分かるが，これは前に求めたAICの値の比較と整合性がある．

ARMAモデルとRandom Forest回帰モデルの比較では，訓練区間では大きな差をつけてRandom Forest回帰のMSEが小さい値となっている．またテスト区間でもRandom Forest回帰が小さいが，差は約７％である．精度的にはRandom Forest回帰が有利であるが，統計モデルを当てはめたARMAモデルでは，予測中央値の他，信頼区間（例えば95%信頼区間）の情報も得られるという点を考慮すべきと思われる．用途に応じて使い分けたい．

（ARMAとRandom Forest回帰のプロット図を見比べてみると，ラインの様子が大きく異なることに改めて驚かされました．また，どのモデルも十分なパラメータ調整ができていない可能性がある点にご注意ください．）


## 参考文献 (web site)
- StatsModels Documentation
http://statsmodels.sourceforge.net/stable/index.html
- Scikit-learn Documentation
http://scikit-learn.org/stable/
- 現場ですぐ使える時系列データ分析（横内大介 (著), 青木義充 (著)）
http://www.amazon.co.jp/ebook/dp/B00KNRL068/
- Python Machine Learning - O'reilly media (By Sebastian Raschka)
http://shop.oreilly.com/product/9781783555130.do