Pythonで学ぶ統計学(13)：予測のための時系列データの統計モデル化（Holt-Wintersモデル）（その２）

前回のブログで時系列データの予測でよく用いられるHolt-Winters法について解説した。
ここでは、pythonを用いて、日本の企業の賃金指数の月次データにHolt-Winters法を適用してみる。
データは、「毎月勤労統計調査 / 毎月勤労統計調査　全国調査 / 長期時系列表月次・年次・年度次・半期・四半期」で、e-statからダウンロードした。
以下は「原系列の図示」「時系列の成分分解」「単純指数平滑法」「二重指数平滑法」「Holt-Winters法」「Holt-Winters法による予測」のコードである。

import pandas as pd
import matplotlib.pyplot as plt
# データの読み込み
wage = pd.read_pickle('wage_index.pkl')
# 系列のプロット (Fig.1)
wage.plot(title='Fig.1 wage series in Japan').autoscale(axis='x',tight=True) 
# トレンド・季節変動・その他に分解 (Fig.2)
from statsmodels.tsa.seasonal import seasonal_decompose 
decompose_result = seasonal_decompose(wage,model='multiplicative')
decompose_result.plot();
span=12; alpha=2/(span+1); wage.index.freq='MS'
# 単純指数平滑法 (Fig.3)
from statsmodels.tsa.holtwinters import SimpleExpSmoothing 
wage['SES12'] = SimpleExpSmoothing(wage).fit(smoothing_level=alpha,optimized=False).fittedvalues.shift(-1)
wage[['w_index','SES12']].plot(title='Fig.3 Statsmodels Holt Winters Exponential Smoothing')
from statsmodels.tsa.holtwinters import ExponentialSmoothing
# 二重指数平滑法 (Fig.4)
wage['DES12'] = ExponentialSmoothing(wage['w_index'],trend='mul').fit().fittedvalues.shift(-1) 
wage[['w_index','SES12','DES12']].plot(title='Fig.4 Statsmodels Holt Winters Exponential Smoothing')
# 三重指数平滑法 (Fig.5)
wage['TES12'] = ExponentialSmoothing(wage['w_index'],trend='add',seasonal='add',seasonal_periods=12).fit().fittedvalues 
wage[['w_index','SES12','DES12','TES12']].plot(title='Fig.5 Statsmodels Holt Winters Exponential Smoothing')
# 直近のみ (Fig.6)
wage[['w_index','SES12','DES12','TES12']].iloc[-36:].plot(title='Fig.6 Statsmodels Holt Winters Exponential Smoothing') 
# Holt-Winters 法による予測
wage = pd.read_pickle('wage_index.pkl'); wage.index.freq = 'MS' # 検証のため、再度データを読み込み
train_wage=wage[:529]; test_wage=wage[529:] # 訓練データとテストデータに分離
fitted_model = ExponentialSmoothing(train_wage['w_index'],trend='add',seasonal='add',seasonal_periods=12).fit() # 訓練
test_predictions = fitted_model.forecast(82).rename('HW Test Forecast') # テスト期間の予測値
train_wage['w_index'].plot(legend=True,label='TRAIN')
test_wage['w_index'].plot(legend=True,label='TEST')
test_predictions.plot(legend=True,label='PREDICTION',figsize=(12,8))
plt.title('Fig.7 Train, Test and Predicted using Holt Winters'); # 実績値と予測値のプロット (Fig.7)