[[時系列データ分析]]
*モデル選択の実際 [#c66db483]
[[Kleiber&Zeileis[2008]>http://www.amazon.co.jp/dp/0387773169]]には、一気に24通りのSARIMAを実行して結果を出すスクリプトが例として載っています(pp.161-162)。このスクリプトはGPL(ver.2)のもとで配布されているAERパッケージに含まれていますから、スクリプトを取り上げて独自の説明をするだけなら著作権侵害にもならないし、本の売れ行きを妨げることにもならないでしょう。この本は計量経済学全体に関するRによって統一された初めての教科書ですし、売れて欲しいと思います。これが売れなければ日本語の教科書なんか出ませんからね。
AERのインストールが終わったとしましょう(2008年以降、並河が授業で配布したR-2.8.0にはインストールしてあります)。スクリプトは
(インストールした場所)\R-2.8.0\library\AER\demo\Ch-TimeSeries.R
にあります。
###################################################
### chunk number 25: arima-setup
###################################################
nd_pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2,
sar = 0:1, sdiff = 1, sma = 0:1)
nd_aic <- rep(0, nrow(nd_pars))
for(i in seq(along = nd_aic)) nd_aic[i] <- AIC(arima(nd,
unlist(nd_pars[i, 1:3]), unlist(nd_pars[i, 4:6])),
k = log(length(nd)))
nd_pars[which.min(nd_aic),]
これが該当部分ですが、これだけでは動きません。(AERとtseriesがインストールしてあれば)単体で動くように補ったものがこれです。
library("AER")
library("tseries")
data("UKNonDurables")
nd <- window(log(UKNonDurables), end = c(1970, 4))
nd_pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2,
sar = 0:1, sdiff = 1, sma = 0:1)
nd_aic <- rep(0, nrow(nd_pars))
for(i in seq(along = nd_aic)) nd_aic[i] <- AIC(arima(nd,
unlist(nd_pars[i, 1:3]), unlist(nd_pars[i, 4:6])),
k = log(length(nd)))
nd_pars[which.min(nd_aic),]
UKNonDurablesはAERに収録されているイギリス非耐久消費財消費額の長期時系列(四半期)データ。AERといっしょに読み込まれるのでdata()で指定するだけで使えます。
先に書いておきますが、UKNonDurablesには強い情報トレンドがあります。見るからに右肩上がりなのです。だから「本当はARIMAでなくてARMAじゃないのか」という可能性は考えないことにします。
先に書いておきますが、UKNonDurablesには強い上方トレンドがあります。見るからに右肩上がりなのです。だから「本当はARIMAでなくてARMAじゃないのか」という可能性は考えないことにします。
では、細かく注釈を入れていきます。
library("AER")
library("tseries")
AERとtseriesをロードします。必要なパッケージも、インストールしてあれば自動でロードします。ただしKleiber&Zeileisの例題群はきわめて多数のパッケージを必要とするためか、AERが自動でロードするパッケージはわずかで、全部を動かしてみるとパッケージ不足で動かない箇所がいくつもあります。
data("UKNonDurables")
複数のデータセットを切り替えて使うときなどに便利な命令。データセットUKNonDurablesの中にある変数名やデータが扱えるようになります。
nd <- window(log(UKNonDurables), end = c(1970, 4))
windowは時系列データの一部を切り出す命令。伸び率が一定だという仮説で、自然対数を取っています。伸び率一定だとなぜ自然対数なのか、などと話し出すとたぶん15分は講義が止まります。
岡部先生の名前も出ているので早苗雅史先生の[[数学玉手箱>http://izumi-math.jp/sanae/MathTopic/MathTopic.htm]]から[[自然対数の図形的意味>http://www.nikonet.or.jp/spring/sanae/MathTopic/e/e.htm]]をご紹介しましょう。ずーっと成長率が変わらない関数って指数関数なんです。そして自然対数を取ると、[[指数関数は一次関数に変換されます>http://ja.wikipedia.org/wiki/%E8%87%AA%E7%84%B6%E5%AF%BE%E6%95%B0]]。
nd_pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2,
sar = 0:1, sdiff = 1, sma = 0:1)
helpファイルによるとexpand.gridは「Create a Data Frame from All Combinations of Factors」。つまりarは0か1か2、diffは1、maは0か1か2…と考えられる3*1*3*2*1*2=36通りのデータフレームを作ります。先走りますが、この時点で
nd_pars
とすると、その内容はこう表示されます。
ar diff ma sar sdiff sma
1 0 1 0 0 1 0
2 1 1 0 0 1 0
3 2 1 0 0 1 0
4 0 1 1 0 1 0
5 1 1 1 0 1 0
6 2 1 1 0 1 0
7 0 1 2 0 1 0
8 1 1 2 0 1 0
9 2 1 2 0 1 0
10 0 1 0 1 1 0
11 1 1 0 1 1 0
12 2 1 0 1 1 0
13 0 1 1 1 1 0
14 1 1 1 1 1 0
15 2 1 1 1 1 0
16 0 1 2 1 1 0
17 1 1 2 1 1 0
18 2 1 2 1 1 0
19 0 1 0 0 1 1
20 1 1 0 0 1 1
21 2 1 0 0 1 1
22 0 1 1 0 1 1
23 1 1 1 0 1 1
24 2 1 1 0 1 1
25 0 1 2 0 1 1
26 1 1 2 0 1 1
27 2 1 2 0 1 1
28 0 1 0 1 1 1
29 1 1 0 1 1 1
30 2 1 0 1 1 1
31 0 1 1 1 1 1
32 1 1 1 1 1 1
33 2 1 1 1 1 1
34 0 1 2 1 1 1
35 1 1 2 1 1 1
36 2 1 2 1 1 1
マニュアルを読んでもわからないことが、こうして中身を出すと直感でわかることがあります。だから手元に環境があることが大事なんです。大学にばかり高度な環境があって、教育が終わったら素手で放り出す、というのは無情なようですが、ソフトを正常流通(利用許諾契約)させるには仕方がありません。Rはその壁を破るソフトなのです。
nd_aic <- rep(0, nrow(nd_pars))
repはrepeat。nrowはnumber of rowsで、行列の行数(いまは36)だけ0を並べたベクトルnd_aicを作れ、という意味。
for(i in seq(along = nd_aic)) nd_aic[i] <- AIC(arima(nd,
unlist(nd_pars[i, 1:3]), unlist(nd_pars[i, 4:6])),
k = log(length(nd)))
help(for)では結果が出てきません。help("for")にするとControl Flow(流れの制御)というhelpページが開きます。
seq(along=nd_aic)というのはnd_aicの要素をひとつずつ処理しろ、という意味のようです。結果的にその数だけ繰り返すことになります。36という数そのものではないところがミソ。単に「36回繰り返せ」でも同じ結果が得られますが、その場合は
for(i in 36)
ではなく
for(i in 1:36)
と書きます。
さて。いよいよAIC関数が出てきました。ndは最初のほうで1970年以後のデータを切り落とした、分析対象のデータでしたね。これをarima関数に食わせて、そこからAICを出そうというようですが…
arima()のマニュアルを読み返しましたが、orderとseasonalを飛ばしていきなり数字の並びを書けば(p,d,q)や(P,D,Q)だと解釈される、とは書いてありません。しかし動いていますから、裏仕様ですね。nd_parsからi行目の1〜3列目、4〜6列目を読んできて、それぞれ(p,d,q)、(P,D,Q)と読ませています。unlistはベクトルや行列の一部を切り出しますが、切り出した後も条件が許せばリストやベクトルになります。
さて。''じつはこのスクリプトではAICは計算できない''といったら驚きますか。驚きませんかそうですか。いや、そろそろ上から読んできて眠くなったかなと思ったので。
最後に「k = log(length(nd))」がついていますが、これはarimaでなくAIC関数の引数です。直前にカンマがあって区切られているのが分かりますね。kを指定しなければ、k=2になってAICが計算されます。こんなkの指定をすると、[[ベイズ情報量基準(BIC)>http://ja.wikipedia.org/wiki/%E3%83%99%E3%82%A4%E3%82%BA%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96]]が代わりに計算されます。
kの指定を削ればAICで大きさを比較することになります。やってみましたが、''AICでもBICでも結局一位は同じ組み合わせでした''。まあ世間ってそんなものですよね。
さて、モデル選択もも大詰めです。
nd_pars[which.min(nd_aic),]
選択はあっけなく1行です。nd_parsは6列ありましたが、ここでは行だけ指定しているので、その列にある6つの数字が表示されるはずです。which.minは読んで字の如く、そのベクトルが最小になる最初の要素番号。
さて、これで解説は終わりました。でもひとつ心残りがあります。sdiffが1に固定で、「季節調整をかけない場合」を比較の対象に入れていませんよね。
まあd=0はありえないとして、pとqであと9通り。ところがexpand.gridの説明に「前のほうに書いた変数から順番に値を回してゆく」とあります。上に挙げた例でもar、ma、sarの順でしょ。じゃ、後ろのほうにダミーくっつけて行列をでっかくして、後ろ半分の先頭9個だけ使って季節調整抜きの分析をやって、最後の一発比較に持っていけばいいんではないか、と思ったわけです。
下にリストを掲げますから、どこが違うか読んでみてください。ちなみに結果は変わらず。季節調整をするとぐぐんと情報量が上がるのがわかりました。
library("AER")
library("tseries")
data("UKNonDurables")
nd <- window(log(UKNonDurables), end = c(1970, 4))
nd_pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2,
sar = 0:1, sdiff = 1, sma = 0:1)
nd_aic <- rep(0, 45)
for(i in 1:36) nd_aic[i] <- AIC(arima(nd,
unlist(nd_pars[i, 1:3]), unlist(nd_pars[i, 4:6])),
k = log(length(nd)))
for(i in 37:45) nd_aic[i] <- AIC(arima(nd,
unlist(nd_pars[i-36, 1:3])),
k = log(length(nd)))
nd_pars[which.min(nd_aic),]
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
