R {seasonal} パッケージで X13-ARIMA-SEATSを使う

X13-ARIMA-SEATS とは

米国国勢調査局 (Census Bureau)で開発している時系列処理/季節調整プログラム。日本の統計データも前身バージョンの X12-ARIMA を利用して季節調整しているのをよく見かける。X13 ではこれまでの季節調整 (X12-ARIMA) に加え、欧州でよく使われている季節調整法である TRAMO-SEATS も利用できるようになった。

(手法自体はちゃんと理解できていないのだが) X12-ARIMA は RegARIMA + 移動平均, SEATS は Wienerフィルタを使って季節性を検出しているらしい。

季節調整法に関する最近の動向：X-12-ARIMAからX-13ARIMA-SEATSへ

X13-ARIMA-SEATS のインストール

Winows

リンク先の winx13.zip をダウンロードして解凍するだけ。

Download Win X-13 - X-13ARIMA-SEATS Seasonal Adjustment Program - US Census Bureau

Mac

ソースからコンパイルするには g77 が必要。gretl のサイトに Mac のバイナリへのリンクがあったので、そっち使う。

http://sourceforge.net/projects/gretl/files/x13as/x13as-x86_64-darwin10.dmg/download

まずはふつうに実行

X13-ARIMA-SEATS はかなり硬派なツールで、入力は独自のフォーマットの .spc ファイルとして渡す必要がある。

Windows

ドキュメント参照 (自分は環境もってないので、、)

Win X-13 Documentation - X-13ARIMA-SEATS Seasonal Adjustment Program - US Census Bureau

Mac

下の tar に入っている testairline.spc ファイルで試す。中身は R の AirPassengers データと同じもの。

https://www.census.gov/ts/x13as/unix/x13asall.tar.gz

実行するには、シェルから以下のコマンドを実行。サンプルファイルを指定する際、拡張子の .spc はいらない。

/Applications/x13as/bin/x13as -i /Users/sin/Downloads/x13asall/testairline

#  X-13ARIMA-SEATS Seasonal Adjustment Program
#  Version Number 1.1 Build 9 PSP = 24
#  Execution began Sat Nov  8 22:35:55 2014
# 
#   Reading input spec file from /Users/sin/Downloads/x13asall/testairline.spc
#   Storing any program output into /Users/sin/Downloads/x13asall/testairline.out
#   Storing any program error messages into /Users/sin/Downloads/x13asall/testairline.err
#  
#   
#  WARNING: At least one visually significant trading day peak has been
#           found in one or more of the estimated spectra.
# 
#  Execution complete for /Users/sin/Downloads/x13asall/testairline.spc at Sat Nov  8 22:35:55 2014

実行結果は同名の .out ファイルに出力される。一部抜粋するが、モデルのフィッティング結果 / 予測値なんかが入っている。

less /Users/sin/Downloads/x13asall/testairline.out

# ...
# 
# ARIMA Model:  (0 1 0)(0 1 1)
#    Nonseasonal differences: 1
#    Seasonal differences:    1
#
# ...
# 
#   Confidence intervals with coverage probability (0.95000)
#   On the Original Scale Before Prior Adjustments
#    ---------------------------------------
#        Date      Lower  Forecast     Upper
#    ---------------------------------------
#    1961.Jan     418.96    442.42    467.19
#    1961.Feb     381.63    412.13    445.06
#    1961.Mar     421.68    463.41    509.26
#    1961.Apr     445.18    496.35    553.40
#    1961.May     439.00    495.82    559.99
#    1961.Jun     501.51    572.95    654.57
#    1961.Jul     579.38    668.98    772.43
#    1961.Aug     558.90    651.83    760.22
#    1961.Sep     471.65    555.17    653.48
#    1961.Oct     411.72    488.94    580.65
#    1961.Nov     351.21    420.60    503.71
#    1961.Dec     395.21    477.08    575.92
#    ---------------------------------------

これはいちいち実行 & 確認がめんどくさいぞ、、。ということで R から使えるようにする。

{seasonal} パッケージのインストール

seasonal というパッケージを使うと、 R から X13-ARIMA-SEATS を呼び出して利用することができる。まずはパッケージをインストール。

install.packages('seasonal')
library(seasonal)

# No path to the binary executable of X-13 specified.
#   
# For installation details, consider Section 2 of the package vignette:
#   http://cran.r-project.org/web/packages/seasonal/vignettes/seas.pdf

ロードすると上記とおり X13-ARIMA-SEATS のバイナリのパスが指定されていないよ！という警告がでる。このままでは使えないので、 vignette に従ってパスを指定する。

補足 Windows の場合は WinX13.exe ではなく x13as フォルダ (x13as.exeを含むフォルダ) へのパスを指定することに注意。指定は <ツールのパス>\\winx13\\x13as になるはず。

Sys.setenv(X13_PATH = "/Applications/x13as/bin")
checkX13()
# Congratulations! 'seasonal' should work fine!
#               - the X13_PATH is correctly specified
#               - the binary executable file has been found
#               - a test run has been successful
# 
# seasonal now supports the HTML version of X13, which offers a more 
# accessible output via the out() function. For best user experience, 
# copy the binary executables of x13ashtml to:
# 
#   /Applications/x13as/bin

R から X13-ARIMA-SEATSを実行

AirPassengers データでモデルのフィット + 予測実行してみる。

# モデルのフィット
m <- seas(AirPassengers)

# フィッティング結果は $data プロパティから取得できるが、このプロパティはすでに deprecate されているため使わないほうがよい
head(m$data)

#         final  seasonal seasonaladj    trend irregular adjustfac
# [1,] 122.7133 0.9019909    122.7133 122.9732 0.9978869 0.9126963
# [2,] 124.7657 0.9542175    124.7657 124.3148 1.0036268 0.9457731
# [3,] 125.0734 1.0698054    125.0734 125.5828 0.9959432 1.0553807
# [4,] 127.5286 1.0023238    127.5286 126.7107 1.0064547 1.0115378
# [5,] 127.3584 0.9486743    127.3584 126.9446 1.0032599 0.9500745
# [6,] 126.1729 1.0762915    126.1729 126.2800 0.9991512 1.0699607

# 代わりに、同じデータを $series から作る
df <- data.frame(final = m$series$s11, seasonal = m$series$s10,
                 seasonaladj = m$series$s11, trend = m$series$s12,
                 irregular = m$series$s13, adjustfac = m$series$s16)
head(df)
#      final  seasonal seasonaladj    trend irregular adjustfac
# 2 122.7133 0.9019909    122.7133 122.9732 0.9978869 0.9126963
# 3 124.7657 0.9542175    124.7657 124.3148 1.0036268 0.9457731
# 4 125.0734 1.0698054    125.0734 125.5828 0.9959432 1.0553807
# 5 127.5286 1.0023238    127.5286 126.7107 1.0064547 1.0115378
# 6 127.3584 0.9486743    127.3584 126.9446 1.0032599 0.9500745
# 7 126.1729 1.0762915    126.1729 126.2800 0.9991512 1.0699607

# 予測の実行
series(m, c('forecast.forecasts'))

# specs have been added to the model: forecast
#          forecast  lowerci   upperci
# Jan 1961 444.2964 418.1670  472.0585
# Feb 1961 413.5093 381.4051  448.3158
# Mar 1961 465.5498 422.4688  513.0240
# Apr 1961 497.5508 445.3449  555.8766
# May 1961 498.3558 440.6152  563.6630
# ...
# Aug 1963 791.9620 502.8166 1247.3808
# Sep 1963 666.3004 417.7790 1062.6581
# Oct 1963 581.5833 360.2171  938.9870
# Nov 1963 510.9198 312.7342  834.6995
# Dec 1963 566.3870 342.6866  936.1157

ツールで実行した結果と違うのは .spc ファイルで行われているオプション指定をやってないからだと思う。

forecast の結果と比べてみる

とりあえずデフォルトの forecast::auto.arima で。

library(forecast)
forecast(auto.arima(AirPassengers), 36)
#          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
# Jan 1961       446.7582 431.7435 461.7729 423.7953 469.7211
# Feb 1961       420.7582 402.5878 438.9286 392.9690 448.5474
# Mar 1961       448.7582 427.9043 469.6121 416.8649 480.6515
# Apr 1961       490.7582 467.5287 513.9877 455.2318 526.2846
# May 1961       501.7582 476.3745 527.1419 462.9372 540.5792
# ...
# Aug 1963       695.2746 574.9455 815.6038 511.2471 879.3022
# Sep 1963       597.2746 473.0909 721.4584 407.3519 787.1973
# Oct 1963       550.2746 422.3523 678.1969 354.6344 745.9149
# Nov 1963       479.2746 347.7200 610.8292 278.0792 680.4700
# Dec 1963       521.2746 386.1853 656.3639 314.6734 727.8759