時系列データの季節調整をしてみよう

時系列データには、季節バイアスが入りがちです。年中行事はもちろんのこと、天候の変化も周期的に発生します。ゴールデンウィークのある5月よりも6月の方が行楽客が少ないとしても、行楽客が減少し出したとは言えないでしょう。そもそも月ごとに日数も異なりますし。月次データを見るときは、季節バイアスの影響を考慮する必要があります。これを数字の処理で行なうのが、季節調整です。実際にRで試してみて、どのぐらい調整できるか見てみましょう。単純移動平均とloessアルゴリズムを用います。

1. 精度を確認する手順

実データでは季節バイアスがどれぐらい入っているのか真実は誰も知らないので、季節バイアスが入った擬似データを作成します。二次方程式から作った非季節変動値に、乱数から作った真の季節バイアスを加えて、観測データを作ります。この観測データを処理して計算された季節バイアスを取り出し、真の季節バイアスと就き合わせてみましょう。

2. 擬似データを作成する

1996年1月から60ヶ月の月次データとして、12ヶ月周期の擬似データを作成します。

# 結果が同じになるように乱数シードを固定
set.seed(20150119)
# 12ヶ月5年分で観測数60
n <- 12*5
# 季節バイアスを乱数から生成
s_bias <- runif(12, min=0, max=4)
# 平均0に正規化しておく
s_bias <- s_bias - mean(s_bias)
# 季節バイアスを連結して、5年分の季節変化をつくる
# tsで時系列データ型にしておく
s_chg <- ts(rep(s_bias, 5), start=c(1996, 1), frequency=12)
# 二次方程式から非季節変動値を作る
x <- 1:n
e <- 0 # 誤差はゼロとしておく。入れるときは rnorm(n) などn個のベクトルを指定。
nsv <- ts(0.9 + 0.003*(x-2*n/3)^2 + 0.1*x, start=c(1996, 1), frequency=12)
# 非季節変動値と季節バイアスを足して、原数値を作成する
y <- ts(nsv + s_chg + e, start=c(1996, 1), frequency=12)

3. 単純移動平均

12ヶ月の算術平均をとる単純移動平均は、単純な計算の割りには効果的に季節バイアスを除去します。最初の11ヶ月はデータが無い欠損値になること、非季節変動値に変化があってから効果が見えるのが遅れことが欠点になりますが、目視するには悪く無いです。

# 移動平均
# 最初11ヶ月分の値はでない。
ma <- ts(numeric(n-11), start=c(1996, 12), frequency=12) # グラフを描くために、時系列データ型にしておく
for(i in 12:n){
ma[i-11] <- mean(y[(i-11):i])
}

4. loessアルゴリズム

Rの組み込み関数stl()で使えるloessアルゴリズムは、局所重み付け回帰関数を用いた高機能な季節調整方法です。ここでは詳しい手順は考えずに、Rに計算させましょう。なお、引数で与える変数が時系列データ型でないと動いてくれないので、注意してください。

stl.y <- stl(y, s.window="periodic")

5. グラフを描いて比較してみる

どれぐらいの精度で季節変動を補正できるのか、グラフを描いて比較してみましょう。この例では、かなりの精度で一致します。移動平均も遅れて動くわけですが、十分機能することが分かります。

plot(nsv + e, ylim=c(0, 10), main="季節調整方法の精度", lwd=2, xlab="", ylab="")
lines(y, lty=2, lwd=1)
lines(y - stl.y$time.series[,1], lty=3, col="red", lwd=2)
lines(ma, lty=4, col="blue", lwd=2)
legend("topleft", c("真の非季節変動値", "原数値", "loessアルゴリズム", "12ヶ月移動平均"), col=c("black", "black", "red", "blue"), lty=c(1, 2, 3, 4), lwd=c(2, 1, 2, 2), bg='white', box.col = 'black', bty="n")

ただし、これは非季節変動値が滑らかな二次関数だからで、ゼロにしている誤差が大きくなると当てはまりが悪くなります。季節変動の標準偏差と同じ大きさの標準偏差を持つ誤差項を入れると、トレンドの変化を掴むのに問題があるほどでは無い*1ですが、真の値と計算値の乖離は目で見て分かる程度になります。

6. X-12-ARIMA

最近はX-13-ARIMAが出てきたようですが、官公庁などでは米国商務省が開発した移動平均型季節調整法であるX-12-ARIMAを使っています。これはアプリケーションが配布されていて、Rからも扱うためのx12パッケージもあります。
Windows機の場合、C:/WinX12にX-12-ARIMAを置いたら、Rでx12パッケージをインストールしたあと、以下のように使います。