![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
[[radvance]] *時系列データ [#je78dfb8] [[そろり新左衛門の倍増話>http://www.nnn.co.jp/dainichi/rensai/naniwa/naniwa070421.html]]は後世の出版による創作話だそうです。t=0の第0日に1粒のコメをもらい、毎日前日の2倍のコメをくれと言ったら、30日後には10億粒ちょっと。「探偵ナイトスクープ」2004年12月10日放送分で数えたところでは1升のコメは96320粒だそうですから、11147升。1升はおよそ1.5kg、コメ1俵は60kgが標準ですから、278俵ほどになります。 このように、時間の経過と共に一定間隔で観測したデータを時系列データといいます。 *時系列データと定常性 [#z02834f2] 「一粒の米の倍増話」は数字がどんどん増えていく時系列データです。ある時系列変数の平均や分散が時間によって変化しないとき、その変数は弱い意味で定常的だといいます。偶然のバラツキを表す撹乱項に定常性がないと、変数の変化のうち、偶然のバラツキ部分とモデルどおりに変化した部分を区別することが出来ません。だから「定常的な部分」を見つけることが時系列分析の鍵になります。 そろり新左衛門は米粒を前日の2倍くれ、といいましたが、別のお伽衆が豊臣秀吉の前に進み出て、前日より2粒多くくれ、といったとします。30日経ってもまだ61粒ですから、1グラム足らず。それでも単調に増え続けているので定常的とはいえません。 でも、次の日もらう米粒との差はいつも2粒で、定常的です。こんなふうに差を取ると定常的になる時系列がある、ということを覚えておいて次の話に進みましょう。 *ARMAモデル [#wca42c3f] ARMAモデルは、経済学で注目される前から、音響の伝達関数などに使われて来ました。残響ってみよーんみよーんと循環しながら減衰したりしますよね。今の音は、1秒前の音や2秒前の音から決まってくるわけです。 あるメーカーにとって、例えばたまたま商品が良く売れたイベントがあるとします。その会社が「最近3ヶ月の平均月間売上」を出すと、翌々月までいつもより売上が高いはず。大きなアクシデントの影響がしばらく残る、ということはよくありますよね。だから、「過去の値」とは別に、「過去の値のうち、ランダムに上下した部分」が今の値に影響している、と考えてもよさそうです。 さっき「1秒前の音」という例を出しましたが、経済問題だと、時間の単位は日や週や月や四半期、ひょっとしたら年でしょう。単位の影響を受けないように「1期前」という言い方をすることにします。Xtの値が、p期前までのXt-pとq期前までのランダム要素εt-qに影響を受けているモデルをARMA(p,q)といいます。 pやqの値を大きくすると、運動の細かいところまで説明できますが、重回帰の変数を増やしすぎるのと同じ問題が出てきます。自由度が下がってくるのです。Xt-pとεt-qの両方を式に入れると、片方だけでモデルを組むより少ない次数(pやqの値)で、運動の大きな部分を説明できることがわかっています。 *さて、どういう推定をするのか [#eb1398d2] ある意味、時系列モデルに意味なんかありません。 例えば日本の貨幣需要がARMAモデル(もう少し後で出てくるARIMAのほうがそれらしいですが、まあ例ですから)に従う、と統計パッケージが示したとします。それで…どうしたらいいのでしょう? 貨幣需要は何から影響を受けるのか。貨幣需要は何に影響を与えるのか。どちらもモデルには含まれていません。 ARMAモデルに代表される時系列分析は[[ルーカス批判>http://hnami.or.tv/d/index.php?%B9%BD%C2%A4%CA%D1%B2%BD%A4%CE%B8%A1%C4%EA%A5%D1%A5%C3%A5%B1%A1%BC%A5%B8%B2%F2%C0%E2#w3f45d56]]に代表される、合理的期待形成を強調する考え方によくフィットしていました。経済変数が自立的に運動する法則があるなら、期待の変化によって構造が変わることはないですよね。金融政策・財政政策の効果について懐疑的な研究結果が次々に現れた時期でもありました。 時系列モデルとはそういうものです。それ自身が経済的に意味のある事実を示すというより、事実として否定できない観測結果を得て、その裏側にある意味や(誰かの)選択肢を考えていくものです。 だから、時系列モデルの形に「妥当な仮定」というものはありません。ARMA(2,2)がARMA(1,1)より必ず良い結果を示すかというと、そうではありません。次数が上がるというのは説明変数が増えるということだからです。 ある程度の範囲の次数をしらみつぶしに試し、[[AIC(赤池情報量基準)>http://ja.wikipedia.org/wiki/%E8%B5%A4%E6%B1%A0%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96]]などの情報量をそれぞれ算出して、最も大きなものを取る、というのが一般的な方法です。 *最尤法 [#f450ced8] &ref(ols.png); のαやβを推定するのに、ふつう使われる最小二乗法は正規方程式と呼ばれる公式を使います。「一番よさそうな値をピンポイントで1回だけ計算する」のです。 これとは別の考え方でαやβの値を推定する方法のひとつが最尤法です。データとして(x1,x2,x3,...xn)(y1,y2,y3...yn)があるとしましょう。αやβを決めると、「x1のもとでy1が観測される確率」を計算できます。厳密には密度関数の値であって確率じゃありませんけれど、その話は置いておいて。なんで確定値じゃなくて確率なのかと言うと、撹乱項εがくっついているからです。 「x1のもとでy1が観測される確率」×「x2のもとでy2が観測される確率」… という数字を計算すると、これは(x1,x2,x3,...xn)(y1,y2,y3...yn)に対してαとβがどのくらい「もっともらしい」かを示す数字と解釈できそうです。これを「尤度関数」と呼びます。 このあたりの解説は、もっとちゃんとした統計学の教科書を読んでください。オンラインで手に入る解説では、[[麻生良文先生(慶應義塾大学)>http://www.law.keio.ac.jp/~aso/]]の講義ノート「[[回帰分析（単回帰）理論編>http://www.law.keio.ac.jp/~aso/ecnm/pp/reg.pdf]]」が明快です。 尤度関数を最大にするようなαとβはどうしたら見つかるでしょうか。αとβを変えて色々試してみる、という力技が、コンピュータの計算能力が有り余っている現在では可能です。実際には、掛け算より足し算のほうが扱いやすいので、尤度関数の自然対数を取った対数尤度関数が最大になるαとβを探します。 AICは、「-2×対数尤度関数の値+パラメータ(説明変数+定数項)の数」です。小さいほど(たいていマイナスの値なので、絶対値が大きいほど)説明力が高く、無駄なパラメータがない、バランスの取れたモデルだと判断できます。 *ARIMAとSARIMA [#fb244cf0] 「差を取ると定常的になる時系列がある」と上の方で書きました。高校では「階差数列」という名前で習っているはずです。ネット上のリソースを探すと、[[北海道算数数学教育会高等学校部会研究部 数学のいずみ>http://izumi-math.jp/]]にある [[札幌新川高等学校 中村文則先生の指導案>http://izumi-math.jp/F_Nakamura/kotewaza/kaisa/kaisa.htm]]に行き当たりました。忘れている人は思い出してください。高校(大学入試)では階差の階差(2次階差)を取るのはよほどの難問だと思いますが、階差数列のそのまた階差数列を取る、というのはイメージできますね。 d階(dは0以上の整数)の階差を取ればARMAになる、というモデルのことを[[ARIMA>http://en.wikipedia.org/wiki/ARIMA]]といいます。0階の階差を取る、というのは何もしないことですから、ARMAそのものですね。だからARMA(p,q)という表記法と同様に、ARIMAはARIMA(p,d,q)と表記します。 さらに、このARIMAに季節調整をかけたモデルがSARIMAです。季節調整の方法は、最近になって飛躍的に難しくなりました。難しくなる前の方法について、[[とやま経済月報>http://www.pref.toyama.jp/sections/1015/ecm/]]2005年11月号「[[統計用語のあれこれ（３）−季節調整について−>http://www.pref.toyama.jp/sections/1015/ecm/back/2005nov/shihyo/index.html]]」 がわかりやすく解説しています。月次データについては前後12ヶ月(厳密には12.5ヶ月)の移動平均(moving average=MA)をとって、それをそのまま季節調整値とするのがよく統計の教科書に載っている季節調整方法でした。日本でも中央官庁のまとめる統計ではESP法、MITI法、センサスX11といった、もっと難しい方法が使われていました。現在では移動平均をARIMAに拡張して、自動で情報量基準を使ってモデルの次数を色々試して一番いいものを取る、というX-12 ARIMA法が導入され、[[日本でもなるべくこれに合わせようということになっています>http://www.stat.go.jp/index/seido/7-1.htm]]。 X-12-ARIMAに関する論文・ソフト・マニュアルの所在については、[[別ページにまとめておきました>http://hnami.or.tv/m2/index.php/%E5%AD%A3%E7%AF%80%E8%AA%BF%E6%95%B4]]。 さて、SARIMAはSARIMA(p,d,q)×(P,D,Q)、あるいはSARIMA(p,d,q)×(P,D,Q)nという形式で表されます。このうちDは「1=季節調整をする、0=しない」だと思ってください。2回以上季節調整をする意味はほとんどありません。nは季節調整の周期で、月次調整なら12、四半期調整なら4です。例えばSARIMA(1,0,1)×(0,1,0)12なら、x(t)のかわりにx(t)-x(t-12)がARIMA(1,0,1)に従うというモデルになります。 Rの場合、SARIMAは関数ARIMAのseasonalオプションで行います。xregオプションを使えば、回帰分析とARIMAを組み合わせたregARIMAも実行できますが、構造変化を含んだモデルにダミー変数を入れるなどの限定的な使い方を除いて、機能が簡単な割りに参考になる情報が少ないので、慎重に使ったほうがよさそうです。 *モデル選択の実際 [#c66db483] [[時系列データ分析(2)]]につづく [[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),]
