Rで季節変動のある時系列データを扱ってみる

Rで計量時系列分析シリーズでだいぶ時系列データの話をしてきたわけですが、最近個人的に季節変動のあるデータを扱うケースが増えてきたので、備忘録的にまとめてみようかなと。

一般に、webデータサイエンスの領域で季節変動というと業種や領域にもよるものの、おおむね

• 週次*1
• 月次*2
• 四半期ごと*3
• 年次or12ヶ月ごと*4

辺りが多いと理解してます（もちろん必ずしもこればかりではないので念のため）。ちなみにこの辺の大ざっぱなまとめが「季節調整」のWikipedia項目に書かれているので、そちらもどぞー。

この辺の処理はRだとかなりお手軽にできるんですが、結構Rならではの約束ごとが多くていきなりやろうとすると「何じゃこりゃ？？？」みたいなことになりがちです。ということで、その辺のポイントをざっくりまとめておきました。

必要なRパッケージ

今回は{forecast}だけインストールして展開しておけばOKです。もちろん、他の関数でも試してみたいという人はCRAN Task View: Time Series Analysisあたりに載っているパッケージの関数で試してみると良いでしょう。

ts(){stats}関数で季節変動を加味した時系列データを作る

Rにはdata(UKgas)など季節変動ありの練習データもあるんですが、この手の四半期データは世の中山ほどあるので、試しにちょっと変な50期周期の季節変動を持つデータを自前で作ってやってみましょう。

> x1<-rnorm(300) # ホワイトノイズ
> x2<-rep(c(15,rep(0,49)),6) # スパイク異常値が6回生じる時系列で周期は50期
> x3<-c(rep(0,150),rep(3,150)) # ステップ状の不連続な変化
> x<-x1+x2+x3 # 全部合わせる
> xt<-ts(as.numeric(x),frequency=50) # 必ずnumeric型で与えること！！！
# なお、frequency=50で50期ごとという周期の情報がインプットされる
> plot(xt)

Rで季節変動のある時系列データを扱う際に、何よりも最も重要なのはts(){stats}関数で生データをnumeric型で渡し、さらにfrequency引数を用いて季節変動成分の目安となる周期インデックスを割り振ってやることです。

なお、生データを与える際はnumeric型でなくてもts()の処理は通りますが、季節変動などのヘッダ情報がNULLになったりカラム数の情報などが正しく渡されないらしく、この後に紹介するもろもろの関数に渡そうとするとエラーの原因になります*5。

stl関数で季節変動を分離する

まず、単純に季節変動を分離したければstl(){stats}関数を使います。使い方は例えば金明哲先生の『Rによるデータサイエンス』にも載ってますし、金先生のサイト（Rと時系列(3)）にも詳しく書かれています。

> xt.stl<-stl(xt,s.window="periodic")
> plot(xt.stl)

こんな感じでバッチリ分離できます。上から生データ、季節変動、トレンド、残差（ランダム成分）の順です。何に興味があるか次第ですが、単に季節調整のみしたデータが欲しいのであればxt.stlからトレンド＋残差のデータを抜き出してやればOKです*6。

季節調整済みARIMAモデルを推定してみる

では、時系列モデリングを実際にやってみましょう。単変量分析の代表例ということで、auto.arima(){forecast}関数でやってみます。auto.arima()にはseasonalという引数があり、これをTrueにすることで時系列オブジェクト内の季節変動情報がロードされます。

> xt.arima<-auto.arima(xt,trace=T,stepwise=T,seasonal=T)
# seasonal引数をTrue (T)にすることで、xtオブジェクト内の季節変動情報がロードされる

 ARIMA(2,1,2)(1,0,1)[50] with drift         : 1010.839
 ARIMA(0,1,0) with drift         : 1537.229
 ARIMA(1,1,0)(1,0,0)[50] with drift         : 1156.663
 ARIMA(0,1,1)(0,0,1)[50] with drift         : 1345.398
 ARIMA(2,1,2)(0,0,1)[50] with drift         : 1233.034
 ARIMA(2,1,2)(2,0,1)[50] with drift         : 983.2024
 ARIMA(2,1,2)(2,0,0)[50] with drift         : 1009.339
 ARIMA(2,1,2)(2,0,2)[50] with drift         : 981.3922
 ARIMA(1,1,2)(2,0,2)[50] with drift         : 978.3383
 ARIMA(1,1,1)(2,0,2)[50] with drift         : 983.5157
 ARIMA(1,1,3)(2,0,2)[50] with drift         : 979.1614
 ARIMA(0,1,1)(2,0,2)[50] with drift         : 1003.731
 ARIMA(2,1,3)(2,0,2)[50] with drift         : 981.5574
 ARIMA(1,1,2)(2,0,2)[50]                    : 976.247
 ARIMA(1,1,2)(1,0,2)[50]                    : 1016.437
 ARIMA(1,1,2)(2,0,1)[50]                    : 976.5955
 ARIMA(1,1,2)(1,0,1)[50]                    : 1015.367
 ARIMA(0,1,2)(2,0,2)[50]                    : 1003.083
 ARIMA(2,1,2)(2,0,2)[50]                    : 979.3051
 ARIMA(1,1,1)(2,0,2)[50]                    : 981.6554
 ARIMA(1,1,3)(2,0,2)[50]                    : 977.052
 ARIMA(0,1,1)(2,0,2)[50]                    : 1001.611
 ARIMA(2,1,3)(2,0,2)[50]                    : 979.4122

 Best model: ARIMA(1,1,2)(2,0,2)[50]                    

Warning message:
In auto.arima(xt, trace = T, stepwise = T, seasonal = T) :
  Unable to fit final model using maximum likelihood. AIC value approximated

50期ごとというかなり重たい季節変動なので少し時間がかかる上に、「これベストモデルじゃないかもよ」という警告が出てしまいました（笑）*7。ともあれこれで季節変動＆トレンドありARIMAモデルが推定されました。forecast(){forecast}関数で推定するとこうなります。

> plot(forecast(xt.arima,range=c(50,95),h=150))

ということで、綺麗に季節変動込み込みで150期先（＝3季先）まで予測することができました。50% / 95%信頼区間についても、季節変動に引っ張られることなくトレンドの部分だけに基づいて推定されていることが見て取れますね。

なお、同様のfrequency / season / seasonalなどの季節調整系の引数は例えばVARselect(){vars}, VAR(){vars}関数などにもついていて、全く同じように季節変動を加味してモデリングすることが可能です*8。

ついでなので四半期データでやってみる

50期周期なんて変なサンプルの例だけで済ませるというのもあんまりなので、普通に四半期データでやってみようと思います。まず四半期データを適当に作ってみます。

> y1<-rep(c(3,5,12,4),20) # 四半期周期のデータ
> y2<-20+20*sin(seq(from=-0.5*pi,to=0,length.out=80))
# 何となく上昇トレンドをsineカーブで表してみる
> y<-y1+y2+rnorm(mean=0,sd=3,n=80) # ランダム成分を加える
> yt<-ts(as.numeric(y),frequency=4) # 四半期周期でtsオブジェクト化
> plot(yt)

ちなみに、ytの中身を表示してみるとこうなります。

> yt
        Qtr1      Qtr2      Qtr3      Qtr4
1   3.000000  5.003953 12.015812  4.035571
2   3.063223  5.098757 12.142159  4.193411
3   3.252493  5.319383 12.394053  4.476474
4   3.566613  5.664435 12.769902  4.882971
5   4.003597  6.131734 13.267330  5.410332
6   4.560683  6.718324 13.883192  6.055223
7   5.234348  7.420496 14.613595  6.813566
8   6.020333  8.233812 15.453919  7.680568
9   6.913669  9.153129 16.398854  8.650747
10  7.908708 10.172635 17.442424  9.717968
11  8.999159 11.285884 18.578031 10.875485
12 10.178127 12.485838 19.798497 12.115979
13 11.438160 13.764911 21.096105 13.431608
14 12.771290 15.115016 22.462650 14.814054
15 14.169090 16.527617 23.889494 16.254577
16 15.622722 17.993783 25.367614 17.744068
17 17.122995 19.504245 26.887667 19.273111
18 18.660424 21.049452 28.440042 20.832040
19 20.225289 22.619636 30.014924 22.410997
20 21.807697 24.204870 31.602356 24.000000

ご丁寧に四半期ということでQtr（四半期）ごとにカラムを並べて、生データを表示してくれます*9。で、上の例と同じようにstl()で分解してみるとこうなります。

> yt.stl<-stl(yt,s.window="per")
# "periodic"は"per"と略記できる
> plot(yt.stl)

ランダム部分に冒されていても、きちんと季節変動が分離されているのが分かります*10。これをauto.arima()で季節調整済み時系列モデルを推定してみるとこうなります。

> yt.arima<-auto.arima(yt,trace=T,stepwise=T,seasonal=T)

 ARIMA(2,0,2)(1,1,1)[4] with drift         : 1e+20
 ARIMA(0,0,0)(0,1,0)[4] with drift         : 440.8788
 ARIMA(1,0,0)(1,1,0)[4] with drift         : 428.1794
 ARIMA(0,0,1)(0,1,1)[4] with drift         : 417.4919
 ARIMA(0,0,1)(1,1,1)[4] with drift         : 418.5298
 ARIMA(0,0,1)(0,1,0)[4] with drift         : 439.311
 ARIMA(0,0,1)(0,1,2)[4] with drift         : 418.5601
 ARIMA(0,0,1)(1,1,2)[4] with drift         : 1e+20
 ARIMA(1,0,1)(0,1,1)[4] with drift         : 417.307
 ARIMA(1,0,0)(0,1,1)[4] with drift         : 416.5326
 ARIMA(2,0,1)(0,1,1)[4] with drift         : 417.5163
 ARIMA(1,0,0)(0,1,1)[4]                    : 428.7213
 ARIMA(1,0,0)(1,1,1)[4] with drift         : 417.3965
 ARIMA(1,0,0)(0,1,0)[4] with drift         : 438.828
 ARIMA(1,0,0)(0,1,2)[4] with drift         : 417.2939
 ARIMA(1,0,0)(1,1,2)[4] with drift         : 418.8026
 ARIMA(0,0,0)(0,1,1)[4] with drift         : 417.5796
 ARIMA(2,0,0)(0,1,1)[4] with drift         : 416.362
 ARIMA(3,0,1)(0,1,1)[4] with drift         : 419.917
 ARIMA(2,0,0)(0,1,1)[4]                    : 426.1714
 ARIMA(2,0,0)(1,1,1)[4] with drift         : 418.1974
 ARIMA(2,0,0)(0,1,0)[4] with drift         : 441.0579
 ARIMA(2,0,0)(0,1,2)[4] with drift         : 418.073
 ARIMA(2,0,0)(1,1,2)[4] with drift         : 1e+20
 ARIMA(3,0,0)(0,1,1)[4] with drift         : 417.9189

 Best model: ARIMA(2,0,0)(0,1,1)[4] with drift         

> plot(forecast(yt.arima,range=c(50,95),h=8))

きちんと季節変動を考慮した上で、右肩上がりのトレンドにあるという予測結果が得られています。

もしかしたら紹介する順番は逆の方が良かったかもですが、日次データなどで季節変動スパイクのあるデータを扱う方が（特にwebデータサイエンス業界では）多い*11と思うので、あえてこの順番で紹介してみました。

最後に

季節変動の周期を間違えるととんでもない推定結果になることがあるので要注意。。。周期の値に不安がある場合は、まず最初にstl()で試しに季節変動を分離してみて妥当かどうか見てから決めた方が良いと思います。

*1:エンターテインメント系のwebサービスは土日にDAUが多かったり

*2:ソシャゲの月初スペンドの跳ね上がりなど

*3:世の中季節が変わるごとに購買行動を起こす人は多いし、政府統計も大半はこれで季節調整をかけてる

*4:年度替わりの購買行動など

*5:stl()にかけた時にエラーが出て「単変量しか扱えません」とかunivariateがどうたらとか言われたらこれ

*6:xt.stl$time.series[,2]がトレンド、[,3]が残差

*7:stepwiseの範囲であまり良い収束が得られなかったということ

*8:ただし関数ごとにかなりのクセがあるので要注意

*9:他の典型的な周期でこのサービスがあるかどうかは不明

*10:それなりにランダム成分が大きくて統計的には微妙なラインかもですが

*11:ただし月次スパイクは正確には30日だったり31日だったり28日だったりするので辛いかも