えんじにあのじゆうちょう

勉強したことを中心にアウトプットしていきます。

コードで理解する移動平均型季節調整法

はじめに

時系列データを扱っていると季節性変動を抽出することがよくでてくる割には、実際どうやってるんだ・・・と思うことが多かったです。
と、言うわけで今回は季節性変動をどのように抽出し、Seasonality, Trend, Residualを分解するかを試してみようと思います。

ちなみに今回は、日本銀行月報5月号を参考にしました。

実装してみる

利用するデータについて

今回もAirline passengersのデータを活用します。理由はわかりやすいトレンドと季節性があるので結果を見るのにいいかなってところです。

df = pd.read_csv("airline-passengers.csv")
df.plot()

f:id:marufeuillex:20200418094153p:plain


見ての通り、右肩上がりの線形トレンドとなんか周期的に上がったり下がったりする波が見て取れますね。

ただ、明らかに後ろに行けば行くほど平均ラインからの分散が増えているので対数をとり、分散をならします。

df["log_Passengers"] = np.log(df["Passengers"])
df["log_Passengers"].plot()

f:id:marufeuillex:20200418094459p:plain

いい感じになったので、このデータを使って進めていきます。

移動平均型季節調整

参考とした日本銀行月報5月号を読むと、概ね以下の手順です。

  • 前提:  y_t = S_t + TC_t + I_t の形に分解する ( S_t: Seasonality,  TC_t: Trend,  I_t: Residual)

1.  TC_t \frac{1}{12} \sum_{i=1}^{12} y_{t-i}とし、Seasonality, Residualを1年分のデータで移動平均を取ることで、排除する。
2.  y_t - TC_tから S_t + I_tを求める
3.  S_t + I_tを月ごとに系列を分ける。そのうえで数年分で加重平均を取ることで I_tは均されるので、 S_tが推定できる
4.  y_t - S_tから TC_t + I_tを求める
5.  TC_t + I_tを加重平均することで、 I_tを均し、 TC_tを求められる
6. 1度目のときは2に戻って繰り返す。2度めなら終了

ということで、順を追ってやってみましょう。

1-1.  TC_tの算出
y = df["log_Passengers"]
tc = y.rolling(12).mean()
tc.plot()

f:id:marufeuillex:20200418094821p:plain

なんか十分トレンドっぽいものが抽出できているように見えますね。
個人的には既に満足感があるのですが進めていきましょう。

1-2. S_t + I_tを求める
si = y - tc
si.plot()

f:id:marufeuillex:20200418095008p:plain

おお、Seasonality + Residualと考えるといい感じに見えますね。

1-3.  S_tを推定

1-2.で算出した S_t + I_tを月ごとの系列を分けます。
そのうえで数年分で加重平均を取ることで I_tは均され、 S_tが算出されます。

加重平均と言われてもどのくらいで取ったものかピンとこなかったので、単純平均で求めてみます。

s = np.zeros(12)
for i in range(12):
    s[i] = si.iloc[i::12].mean()

plt.plot(np.tile(s, tc.shape[0]//12))

np.tileは与えた配列的な要素を指定した個数分並べるメソッドです。今回sは12しかないので、tcの数に合わせて繰り返すことにします。

f:id:marufeuillex:20200418095607p:plain

いい感じですね。(それしか言ってない)

1-4.  TC_t + I_tを求める

 y_tから1-3.で求めた S_tを引くことにより TC_t + I_tが求められます。

ti = np.zeros(tc.shape[0])
for i in range(12):
    ti[i::12] = y[i::12] - s[i]
plt.plot(ti)

f:id:marufeuillex:20200418095904p:plain

確かにTrendに残差が乗っているように見えますね。

1-5.  TC_tを求める

1-4.で求めた TC_t + I_tを加重平均することで、 TC_tを求められます。
また、どのくらいで加重平均を取ったものか悩ましいので単純平均にしてみました。

tc2 = pd.Series(ti).rolling(12).mean()
plt.plot(tc2)

f:id:marufeuillex:20200418100229p:plain

おおっ、いい感じ・・・なんでしょうか?なんだか1-1.で求めたものと差がない気がします・・・
確認のため、1-1で求めたものとyを並べて見てみます。

plt.plot(tc, label="tc")
plt.plot(tc2, label="tc2")
plt.plot(y, label="src")
plt.legend()

f:id:marufeuillex:20200418100416p:plain

青が1-1、オレンジ色が1-5で求めたものです。
いいかどうかはわからないですが、たしかに若干異なってますね。そして、どちらにせよyのトレンドを表せているように見えます。
続けていきましょう。

次は手順2に戻って繰り返します。
以下は同じ手順なので、淡々と進めます。

2-2. S_t + I_tを求める
si2 = y - tc2
si2.plot()

f:id:marufeuillex:20200418100804p:plain

2-3.  S_tを推定
s2 = np.zeros(12)
for i in range(12):
    s2[i] = si2.iloc[i::12].mean()

plt.plot(np.tile(s, tc2.shape[0]//12), label="s")
plt.plot(np.tile(s2, tc2.shape[0]//12), label="s2")
plt.plot(y, label="src")
plt.legend()

f:id:marufeuillex:20200418101042p:plain

2-4.  TC_t + I_tを求める
ti2 = np.zeros(tc2.shape[0])
for i in range(12):
    ti2[i::12] = y[i::12] - s2[i]
plt.plot(ti2)

f:id:marufeuillex:20200418103725p:plain

2-5.  TC_tを求める
tc3 = pd.Series(ti2).rolling(12).mean()
plt.plot(tc, label="tc")
plt.plot(tc2, label="tc2")
plt.plot(tc3, label="tc3")
plt.plot(y, label="src")
plt.legend()

f:id:marufeuillex:20200418103813p:plain

結果の確認

では、Trend, Seasonality, Residualを並べて確認してみましょう。

i = tc3 - ti2
fig, ax = plt.subplots(4, 1)
ax[0].plot(y)
ax[0].set_title("src")
ax[1].plot(tc3)
ax[1].set_title("Trend")
ax[2].plot(np.tile(s2, tc3.shape[0]//12))
ax[2].set_title("Seasonal")
ax[3].plot(i)
ax[3].set_title("Residual")

f:id:marufeuillex:20200418101807p:plain

おお、よくできてそうですね。

最後に重ね合わせてみます。

y_hat = tc3 + i
for j in range(12):
    y_hat[j::12] += s2[j]
plt.plot(y_hat, color="blue", label="pred")
plt.plot(df["log_Passengers"], color="red", label="y")
plt.legend()

f:id:marufeuillex:20200418101959p:plain

んー、まぁ今回加重平均使ってなかったり適当にやってるのでこんなもんか、っていう感じですが、手順の確認という意味では概ね大丈夫な気がしますね。

ライブラリでやるとどうなるの?

statsmodelsのライブラリにseasonal decompositionの関数があるので同様に試してみましょう。

import statsmodels.api as sm
res = sm.tsa.seasonal_decompose(y, period=12)
res.plot()

f:id:marufeuillex:20200418102124p:plain

簡単ですね。

一応重ね合わせてみます。

seasonal = res.seasonal.dropna()
trend = res.trend.dropna()
residual = res.resid.dropna()

y_hat = trend + residual
for j in range(12):
    y_hat[(j-6)%12::12] += seasonal[j]
plt.plot(y_hat, color="blue", label="pred")
plt.plot(df["log_Passengers"], color="red", label="y")
plt.legend()

f:id:marufeuillex:20200418102325p:plain

おおっ、完璧に重なってますね。ちゃんとやればこうなるっていうことで勉強になりますね。

終わりに

今回は季節調整法の最も基本的な移動平均型の季節調整についてまとめてみました。
自分で実装することにより、データに対してだいたいこういうことをやるということがわかり、実際のデータに適用するときにもよく理解しながら利用できますね。

あと、今回データは対数をとっているので、最終的に使うときはもとに戻してお使いくださいね。