Edited at

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

More than 1 year has passed since last update.

(追記)

本記事後半部,Random Rorest Regression に関する訂正記事を投稿しました.こちらも参照ください.時系列データ分析の処理でやってはいけないこと(反省を含めて)

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 がサポートされていれば,上記4つのモデルはすべて対応可能となる.但し,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

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


ARMAモデルの次数選択

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

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

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

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) が良さそうとの情報を得た.よってこの2つのモデルをデータにフィットさせてみる.

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) も試してみた.

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の値を確認する.

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年)のモデル推定値を算出する.

# 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クラス)のメソッドで求めることができるような実装となっている.元データの値,推定値をまとめてプロットする.

# 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

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

# 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

階段状の直線でフィットすることが特徴となっている.今回は単変量の時系列データであるが,いくつかの過去のデータ使って現在の値を推定するモデルとしてみた.具体的には,

Volume_current ~ Volume_Lag1 + Volume_Lag2 + Volume_Lag3

というモデルである.訓練データとテストデータを前回と同じ(70, 30)にしたかったが,Lag値を計算する過程でNaNがでるので,初めの部分をdropna()して,(67, 30) のデータ長とした.

まず,データの前処理を行う.

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を使用.

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

最後にMSEを算定する.

# 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回帰が小さいが,差は約7%である.精度的にはRandom Forest回帰が有利であるが,統計モデルを当てはめたARMAモデルでは,予測中央値の他,信頼区間(例えば95%信頼区間)の情報も得られるという点を考慮すべきと思われる.用途に応じて使い分けたい.

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


参考文献 (web site)