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

はじめに

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

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

実装してみる

利用するデータについて

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

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

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

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

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

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

移動平均型季節調整

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

• 前提: の形に分解する (: Seasonality, : Trend, : Residual)

1. をとし、Seasonality, Residualを1年分のデータで移動平均を取ることで、排除する。
2. からを求める
3. を月ごとに系列を分ける。そのうえで数年分で加重平均を取ることでは均されるので、が推定できる
4. からを求める
5. を加重平均することで、を均し、を求められる
6. 1度目のときは2に戻って繰り返す。2度めなら終了

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

1-1. の算出

y = df["log_Passengers"]
tc = y.rolling(12).mean()
tc.plot()

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

1-2.を求める

si = y - tc
si.plot()

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

1-3. を推定

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

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

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の数に合わせて繰り返すことにします。

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

1-4. を求める

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

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

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

1-5. を求める

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

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

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

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

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

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

2-2.を求める

si2 = y - tc2
si2.plot()

2-3. を推定

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()

2-4. を求める

ti2 = np.zeros(tc2.shape[0])
for i in range(12):
    ti2[i::12] = y[i::12] - s2[i]
plt.plot(ti2)

2-5. を求める

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()

結果の確認

では、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")

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

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

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()

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

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

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

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

簡単ですね。

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

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()

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

終わりに

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

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