ログインしてさらにmixiを楽しもう

コメントを投稿して情報交換!
更新通知を受け取って、最新情報をゲット!

Rで経済時系列コミュのAIC(赤池の情報基準量)

  • mixiチェック
  • このエントリーをはてなブックマークに追加
#お約束?のAICです。連載っぽくしようか?と思っていましたが、一回の読みきりにしました。

生真面目にAIC(赤池の情報基準量)について話をしようか?とも少しは思ったのですが、いつものごとく、厳密な扱いに関しては何かの教科書で学んで頂くこととして、ざっくりAICの話をして、まがいなりに?ではあっても、実用上扱えるようになる事を目指す事にします。

ざっくりAICとは何か?と言えば、実際にモデルはどれを選択すれば良いのか?と言う問いに対して、明確な答を出す方法の話です。

複数のモデルがあった場合、どっちのモデルを選んで使えば良いのか?と言えば、当然の話として、より優秀なモデルを使えば良いに決まっている話しなのですが、ではそもそも優秀なモデルとは、どういう事よ?と言うのを旨く説明できるでしょうか?

例えば、身長を予測するとします。この場合、体重、スリーサイズ、性別、年齢、など色々な予測に有用なデータが、ありそうに思えます。

出来るだけ正確に予測しようとした場合、なんとなく、こういったデータの種類が多いほど良さそうな感じもしますよね?

実際に、こういったデータを増やせば増やすほど誤差は減っては来るのですが、意味があって誤差が減ったのか、やみくもに見掛けだけの誤差が減ったのかが、判らないのが困った話なわけです。

では、AICを使う事で、どうしてこの問題が解決するのかと言えば、AICは誤差が減る事も評価するのですが、何種類のデータを使うのかという事も評価に使っている為、誤差の少なさと、使う種類の少なさの、双方の均衡の取れたモデルを選択出来るようになるという事です。

纏めると AICは

1)出来るだけ誤差が少ないほうが良い
2)できるだけ説明に使うデータの種類の少ないほうが良い

と言う相反する問題をバランスさせて評価を決定する方法な訳です。

ここで、以下のような実験をしてみました。データは檜木さんのデータを
使い、対数正規過程だとしてARモデルを当て嵌めてAICを眺めている所です。

> n<-read.csv("nkf.csv")  データを読む
> d.l.n<-diff(log(n[,2])) 対数変換し和分する
> model <-ar(d.l.n)    ARモデルを作る
> model$aic        AICを表示させる
0 1 2 3 4 5 6
2.343794 0.000000 1.990054 3.359000 3.145408 4.704977 6.700395
7 8 9 10 11 12 13
6.790411 7.718905 3.827878 4.186501 6.076982 8.024328 9.185088
14 15 16 17 18 19 20
10.278655 11.880014 11.356852 10.857688 8.334711 9.720635 9.428875
21 22 23 24 25 26 27
11.232674 13.206178 14.517388 15.011147 15.926222 16.733067 18.711458
28 29 30
17.641039 16.714192 18.657771
>

ARの次数とAIC値が表示されているのですが、一番値が小さいのは次数が1の時でAIC値は0です。(AIC値は小さい程良い値です。)

つまりAICは、AR(1)を選択するのがベストだと言っているわけです。

ここで、次数が0の時の値が 2.343794 と出ています。AR(0)モデルという事なのですが、これがランダムウオークモデルです。

計算する上でAR(0),MA(0),ARMA(0,0),ARIMA(0,1,0)モデル等は、ランダムウオークモデルだという事です。

また次数が2の時のAIC値は1.990054ですので、AR(1)よりは劣るがAR(0)(ランダムウオークモデル)よりは、優れている事を示しています。

実際に金融予測モデルを立てる場合には、将来の数値が、上がるのか下がるのかが問題になる訳で、モデル選択的に重要な話は、ランダムウオークモデル(上げ下げが期待値を持たないモデル)よりも、モデルとして優れていて初めて、上げ下げ予測としての価値があると考えられるという事です。

つまりこの問題で言えばAR(3)以上はモデル選択論的に無意味なモデルと考えて良いという事です。

ざっくりとこんな感じでAICを使う為、金融に限らずモデラーはモデルを作ってはAICばかり計っているのが実質の、お仕事という話です。

コメント(23)

ここでやってるモデル類のAICの実装ですが、半ば暫定的に以下で計算することにします。

#一般式とは違う形式だけど、最小二乗の場合は、こうなる筈だと。。
#まぁ、これで間違ってると言う事も無いだろうと思いますが。

AIC = n(log(2πS/n)+1)+2(p+2)

n: サンプル数
S: 残差平方和(二乗誤差を全部足したもの)
p: モデル次数

参考:
http://software.ssri.co.jp/statweb2/tips/tips_10.html

で、取り合えず実装例


--------------------------------------------------------------------
library(matlab)
library(signal)
Fs <- 1000;N <-150; time<-1:N; N2<-2; N3 <-2;
a <-matrix(c(1.4 ,-0.98),nrow=1,ncol=2);
a2<-matrix(c(-0.98, 1.4),nrow=2,ncol=1);
aa<-c(1 , -a);
Nrn<-rnorm(N);
yobs<-matrix(filter(1,aa,1*Nrn),1);
#par(mfrow=c(3,1));#画面を3×1 に分割
#plot(time,yobs,type="l");axis=c(0, N,-20, 20);

# State Space Representation
#
# x(t+1) = F x(t) + G w(t) # cov{w(t) w(t)} = {Q 0}
# y(t) = H x(t) + v(t) # {v(t),v(t)} = {0 R}
F <- matrix(c(1,0,0,1),2,2);
G <- matrix(c(1,0,0,1),2,2);
H <- matrix(c(1,1),1,2);
Q <- matrix(c(1,0,0,1),2,2);
R <- matrix(c(1),1,1);
# ----------------------------------
# Kalman Filter
# ----------------------------------
#
# initial value
x <- matrix(0,N2,1); # x(0|0)
P <- diag(N2); # P(0|0)
ypre <- matrix(numeric(N),1);
Spre <- matrix(numeric(N),1);
K <- matrix(0,N2,1);
Start <- 3;
ypErr2 <- 0;
#
for(t in Start:N){
#
# prediction
x <- F %*% x; # x(t+1|t) <- x(t|t)
P <- F %*% P %*% t(F) + G %*% Q %*% t(G); # P(t+1|t) <- P(t|t)
H<-matrix(c(yobs[,t-1],yobs[,t-2]),1,2);
ypre[t] <- H %*% x;
Spre[t] <- sqrt(H %*% P %*% t(H)+R);
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R); # Kalman Gain
ypErr <- yobs[,t] - ypre[t];
ypErr2 <- ypErr2 + ypErr^2;
x <- x + K * (ypErr); # x(t|t) <- x(t|t-1)
P <- P - K %*% H %*% P; # P(t|t) <- P(t|t-1)
}
#plot(time,yobs,type="l");
#par(new=T)
#plot(time,ypre,type="l",col="red",xlab="",ylab="",axes=F);
#plot(time,Spre,type="l");

AIC <- (N-Start)*(log(2*pi*ypErr2/(N-Start))+1) + 2*(N2 + 2);
試してないんでわかりませんが、引数を受け取った後にabs()するのでは駄目なんでしょうか?

仮にマイナスのノイズを渡されても、中でプラスになるわけで??
>確かペナルティ関数がどうたらというのあったなーと思って。

ええ、正統派のやりかたで、ABS()なんてほうが、インチキくさいわけで。。

まぁけど、常々思っている事だけど、両方必要な事なんですよ。

考えてみれば、ノイズの分散値にマイナスが同定されることだって、数学的に見れば虚数ノイズが載っていると言っているわけだし、じゃ、これどういう意味があるんだろうか?とか考えるのも、無意味な話でもないわけですよね?

まぁけど、隠れモデルの推定変数(回帰係数など)が、分散を持つって事自体、自然に感じます?

ベイジアンなら当然の話ではあるけれど、最初は気持ち悪かな?と。

カルマンの動的推定は、窮屈な限定事項はあるけれど、ベイズ推定の一種でもあるわけです。

相場だ!カルマンだ!と、オモチャを色々並べて遊んでおいて、こういう本質が見えてくるのが面白いなと感じているところ。
>レトロ系シストレの人は「カーブフィッティングに陥らないように」なんてことを考えつつ、利益が良くなるようにパラメータを調整してしまいます。

つか、最適化と適応化の違いは、まさにこの部分なわけで。。

調整するとしたら、過去の相場で利益が出たようにしか調整できないし、そうするとそれが明日の相場には使えないわけでね。。

ここでやってる事は、本来はシストレと言うよりも、シストレの入門部分という事です。

「非定常な時系列データの動的な推定」が出来ないならば、経済学、金融工学、シストレ、全部合わせても当然ペケです。
>今のところ、あまり賢い探索になってない予感はします。
>実行速度のわりにそれほどいい結果が得られていないという。

 関数kalmanのコードの日付が一日進んでいるのが原因では?

 『ypreの日付進み』等を修正しましたら、
 非常に『賢い探索』になっているというのが、私の実感です。

 温泉先生の言葉を借りれば
> トレンド推定はしているものの、predictionの中では、
> 実は1期前の推定値がそのまま出てきているだけ。 

 そもそも、OctaveからRコード変換したときの日付進みミスを
 いままで気づかずに使用してきたことに原因があり、
 私のコード変換ミスを痛感します。大変申し訳ありませんでした。
>使えると信じている人はそっとしておいたほうがいいのかな、とも思うんですが。

証券業界って、つくづく色々巧みだよなって思うわけです。

計算なんかしなくても、相場がダイナミックな事くらい誰でも判っていること。

で、そういうの相手に14日のRSIが30以下なら買え!って、スタティックな事を言って意味ある筈も無いわけで。。

競馬はJRAが50%のテラ銭とるけど、相場はGSが新入社員にも1千万以上のボーナス出したりするわけで、さてどっちが、トータルだと、きついんですかね?。。(ぼそ
レヴィさん
>私のコード変換ミスを痛感します。大変申し訳ありませんでした。

いや、気にすることは、全くないです。

なんだかわからないから、勉強するのだから、途中で間違うのは
寧ろ当然だと思いますよ。

檜木 さん
>オシレータがどうとか、ノイズにフラれるものを信じてしまう。

今回、トレンド項を入れたのは、良かったか悪かったかどうかな?とは今でも思っている事なのですが、今回のような雑なやりたかではあっても、トレンドのシステムノイズは推定できると思います。

ボラの数倍あるようなトレンドノイズが同定されたら、推定誤差の問題では無くなるわけです。

まぁ本当は、以下の大先生のモデルなんかでやるほうが良いのでしょうけどね。

http://www2.ics.hawaii.edu/~gersch/

#欲しい論文などがあれば、かっさらっておいてください。
#昨年ご他界なさったらしく、いつクローズされるかわかりません。
>いや、気にすることは、全くないです。
>なんだかわからないから、勉強するのだから、途中で間違うのは
寧ろ当然だと思いますよ。

なにを偉そうに、オオボケを言ってるんだろう。。と、ちょいがっくり。(自分

そもそも、これはレヴィさんの、間違いでもなんでもないんです。

回帰係数は、予測の話、トレンドは現在値推定の話です。

たまたま運良く、トレンドの現在値推定を、現在ではなく昨日のものにしてあったのが幸いしているわけで、トータルでは予測になっているわけです。

Yt = Yt-1 + a1 x (Yt-1 - Yt-2) + εt

とすれば、これは和分したAR(1)です。今回precompと言ってるものは

Yt = Tt-1 + a1 x (Yt-1 - Yt-2) + εt

となっている為、トータルでは予測になっているわけですが

Yt = Tt + a1 x (Yt-1 - Yt-2) + εt

としてしまうと、これはトレンド推定の不足分を予測で補う現在値推定になってしまう為、予測モデルではなくなるわけです。

トレンド項は、カルマン推定されFを通過した時点で、決定しているわけで、回帰は次の予測のHに出会って時系列の値を貰って初めて値をもつわけです。

と言う事で、Cで書いたソースは現在値の推定だけをやるなら、予測風に一期遅らせるのではなく、Fを通過する時間(つまり推定と同一時間)を、値として持たせているわけです。

まぁ予測風にした問題は、和分したのと比較して、予測モデルとして、どうメリットがあるか?と言うのが不明なところなのですが、トレンドのシステムノイズを、不細工ながらも推定可能になっているところが、和分したものよりも、メリットがあるのでは?と考えていると言う事です。

いやこんなもの要らないやって人は、ARの和分モデル風に記述すれば、そのまま非定常な時系列が扱えるわけです。

と言う事で、私がRでトレンドモデルを作った時の時系列的な時間の扱いを間違ったと言うのが真相で、これが運良く予測モデルになっていたと言う事です。

申し訳ありませんでした。
まぁこれでいいか?と思った、もうひとつの理由は次数なんです。

ARIMA(1,1,1)モデルであれば次数は3ですがTARMA(1,1,1)でも次数は3なわけです。

本来は、両方作ってAICでも比較すれば良い話なわけですが、少なく見ても和分よりは器用そうだし、それでいて次数が同じなら、儲けかも?と言う程度で考えています。

ここのトレンドモデルには極がある為、Q=R=0に出来ませんが、こうすれば意味は和分なわけです。

まぁQ,R付きなのに和分と同じ次数でカウントしていいのか?と言うのは、ちょいインチキくさくも感じますが、、、
温泉先生
いつもお世話になり有難うございます。
復習と訂正を兼ねての確認質問をさせて下さい。

 ?日付進みに関する(prediction)と(filtering)部分の誤字訂正(2箇所)は
   以下で正しいのでしょうか。
   その他、訂正すべき箇所があればコード訂正指示をお願いします。
   よろしくお願いします。

F <- matrix(c(1,0,0,1),2,2);
G <- matrix(c(1,0,0,1),2,2);
H <- matrix(c(1,1),1,2);
Q <- diag(2);
R <- matrix(0, 1,1);
# ----------------------------------
# Kalman Filter
# ----------------------------------
# initial value
#
x <- matrix(0,N2,1);
P <- diag(N2);
ypre <- matrix(numeric(N),1);
Spre <- matrix(numeric(N),1);
K <- matrix(0,N2,1);
Start <- 3;
End<-N+1;
ypErr2 <- 0;
#
for(t in Start:End){
#
# prediction
#
x <- F %*% x;
P <- F %*% P %*% t(F) + G %*% Q %*% t(G);
H<-matrix(c(yobs[,t-1],yobs[,t-2]),1,2);
#ypre[t] <- H %*% x;  誤
ypre[t-1] <- H %*% x;#正
Spre[t] <- sqrt(H %*% P %*% t(H)+R);
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t(H)+R);
#ypErr <- yobs[,t] - ypre[t];   誤
ypErr <- yobs[,t-1] - ypre[t-1];#正
ypErr2 <- ypErr2 + ypErr^2;
x <- x + K * (ypErr);
P <- P - K %*% H %*% P;
}
n<-End-Start;
AIC <- n*(log(2*pi*ypErr2/n)+1) + 2*(N2 + 2);
AIC


>ここのトレンドモデルには極がある為、Q=R=0に出来ません

 ?『極がないトレンドモデル』とは
   どのようなRコード記述になるのでしょうか?

一方的な質門ばかりで恐縮しますが、
以上よろしくお願いします。
>?日付進みに関する(prediction)と(filtering)部分の誤字訂正(2箇所)は

いえ、訂正が必要ないと思っているんですけど?

ループの中でt期の予測を作り、t期の実値と比較(答えあわせ)し誤差を勘案し、次のループで使う回帰係数を作り

と言うのが一周目で

二周目からは、貰った回帰係数で予測し、、

と言う手順を繰り返しているので、訂正が必要ない筈。

もしも訂正してしまうと

Yt = Yt + Yt-1 + εt

と言うモデルになるだろう?と思いますけど?

>?『極がないトレンドモデル』とは

と言うか、どうして極を持つかなんですが、私のCのソースで言えば

// filterling
K = P / (P + R);

とカルマンゲインの計算そのものが、P=R=0では、計算できないわけです。
と言うかカルマンフィルタ自体が、誤差の正規性がある事が前提で
計算されるため、計算できなくなる為です。




温泉先生
早々の回答有難うございます。
理解しやすい回答に感謝します。


(上記2の、取り合えず実装例)に 
17の?日付進みに関する誤字訂正(2箇所)をしたRコードで
AICを計算しました。

>AIC
[1] -130.8305

ここで、次の質問で恐縮ですが、
マイナスのAIC計算結果に戸惑い、かつ
どのように評価すればよいのかわかりません。

上記の質問に関する回答よろしくお願いします。
>マイナスのAIC計算結果に戸惑い、かつ
>どのように評価すればよいのかわかりません。

まず、評価ですが、小さい程よいと判断するでOKです。

問題なのは、AICがマイナスなところ。

AIC = -2 ln(L) +2K

が一般式なので、AICがマイナスになると言う事は
少なくともln(L)がプラスである必要がある。
つまりL>1である必要があることになる。
ここでLは尤度なので、確率を逆さまにみると1を
超えている事になる。

と言う事になれば、戸惑うわけですよね?

さいころを投げて、1の目が出る確率は1/6で
これを逆さまにみると、面が1/(1/6)で6面ある
事になる。
もしも面が6面以上あるのが正常と言う事になれば
そもそも推定された確率(1/6)そのものがおかしい。

だからこのAICの計算は、おかしいんじゃないか?と。。

と、そういう言う話でいいですか?

温泉先生
判りやすい説明有難うございました。

>まず、評価ですが、小さい程よいと判断するでOKです。

 今後、マイナス数値の大小は、誤差僅少の程度(良適合の度合い)
 であると評価するよう考えを改めます。有難うございました。
いや、自分でも計算して、話のツジツマが合わないので困っていて
調べていたんです。

尤度を確率そのものだと考えると、0〜1の間しか値が取れないわけで
ふに落ちないんですよ。

ところがRのARなんかを見ても最小AICが0になるように調整されている
んです。

まぁ、モデルの良さそのものは、残差平方和が小さいほど良いだろうけど
この時に同じことなら、次数が小さいほうが良いわけで、結局さまざまな
基準量の論争も、回帰であれば次数と残差平方和を、どうコンペアすれば
良いか?という話になるわけなのだけど、動的予測モデルの場合は最終的
な判断が、結局向きがあたるかどうか、つまり儲かるかどうかと言う事に
なる為、相場に限らずμ/σが問題になるわけです。(天気で言えば晴れか
雨が当たれば良い)

但し、次数を上げると、動的モデルの場合σそのものが、どこからか悪化
するので、μ/σも、悪化を始めるのですが、これは回帰で言えば係数の
動的追従そのものが次数を上げると悪化する為と思われます。

勿論、σが上昇すれば、残差平方和が悪化しているわけで、これはどの
基準量を採用したところで、悪化している事になるわけで、次数との
コンペアの問題ではなくなるわけです。

こんな事を平気で言えば、怒られるんだろうけど、そもそも本質的に
モデル次数に絶対的な意味があれば、AIC以降に発表された他の基準量が
幅を利かせている筈なのですが、色々な主張はあっても、そうはなって
いないわけです。

この事は少なく言って見れば、誰も本質的な次数と尤度の関係そのものは
実際には問題にしていないと言う事です。

もしも本質的な意味があるのであれば、そういう新しい基準量に移っているか
やっぱり本当はAICだけが正しかったかと言う事になるのですが、実際には
複数の基準量が記載されたモデル資料なんかを見るのもある話なので、つまり
まぁ一種のモデルの効率が表示してあるだけと見ておくべきだと考えています。

また、こないだのQとR付きでも、これがモデル次数に勘定されないが本当に
正しいのであれば、これらはモデルの問題ではなく、本質的にもっていた
ノイズだと考える事になるわけですが、これだって良く考えてみると
モデルを作らなければ、少なくともQのモデルノイズは無いわけで、やはり
おかしな話ですが、動的モデルだから次数が上がっているってのも逆に
見たこともないんです。

つまりモデル選択論を、厳密に言えば、結論が出ていない分野でもあり
相対的な評価のペナルティの課し方が違うだけのことで、その絶対値が
何か実量的な意味を持っているとは考え難いと判断しています。

RのARのAICなんかも、結局はそういう範囲であれば最小を0にすると
言うのは、可読性が上がる方法だと言う事だろうと判断しています。

まぁ、残差平方和そのものが小さくなるほうが良く、それと次数の関係の
ペナルティ割合が、AICになっていれば、使うぶんには問題ないわけだし
ということです。
以下は私見だし、こんな事を言う人が他にいないので聞き流して
欲しいんですが、

例えばカルマンフィルタを使うのあれば、誤差は正規性を仮定してる
わけで、まずはこの時点でモデル選択数1が入り、QやRを使うので
あれば、ここにも1が入るので、選択数は最小でも3になる。

って誰も言わないし、私もそうだと思わないんですが、じゃこれらを
仮定しなかったら、話そのものが成立しないわけだし、変数が増える
ほどQの数が増えるわけで、これも無視してるわけで、やっぱおかしい
と思うわけです。

では、どうしてそうしないのか?と言われれば、

1)みんなそうしているから
2)相対的な評価法だと思っているから

です。

ところがAICでもKL情報量からの直伝みたいな話になっているしBICなんかは
ベイズ推定上こうなる筈だと書いてある。

で、両方ともQやRはカウントしない。

こりゃ少なくとも自分としては相対的な評価指数みたいにしか使ってはいけない
ものだな。。と思っていると言う事。

まぁ、どの基準量も根拠がきちんとあると主張していますけどね。。(ぼそ

ログインすると、残り6件のコメントが見れるよ

mixiユーザー
ログインしてコメントしよう!

Rで経済時系列 更新情報

Rで経済時系列のメンバーはこんなコミュニティにも参加しています

星印の数は、共通して参加しているメンバーが多いほど増えます。

人気コミュニティランキング