★新サイト完成しました!
3秒後に自動的に移動します
変わらない方は こちらからどうぞ
http://logics-of-blue.com/%E6%99%82%E7%B3%BB%E5%88%97%E8%A7%A3%E6%9E%90_%E7%90%86%E8%AB%96%E7%B7%A8/
ここからは、予測理論 の本命、時系列解析の紹介に移ります。
時系 列解析は、回帰分析とは違ってあまり知らない人も多いと思うので、ざっと解説を載せておきます。
それ と、便利なパッケージ forecast の紹介も。
回帰分析だと、なにかの要因から結果を予測します。一方時系列解析は、過去のデータから未来を予測するという手法です。
たとえば、一年前に魚が豊漁だったら、今年も豊漁になりやすいとか、そういう関係性があるのかというところを調べて予測に活用します。もちろん、回帰分析
のように別の要因の影響を取りこんで予測モデルを作ることも可能ですが、とりあえずここでは一変量のデータのみを扱います。ようするに、過去の自分のデー
タしか使わないということです。
まず、時系列データにおいて、回帰分析は余りお勧めできません。
なぜか。
例をあげます。たとえば、北海道に牛の赤ん坊が生まれて、また沖縄にも同時刻に人間の子供が生まれたと。このとき、牛の子供と人間の赤ん坊の体重を回帰
させてやると、おそらくきれいな回帰直線がひけます。両方とも成長しているからですね。
でも。そんな回帰分析には何の価値もないでしょう。北海道の牛がダイエットしたから沖縄の少年の体重が軽くなるなんてことはあり得ないわけです。でも、
回帰直線は引っ張れてします。ようするに、誤った予測をしやすくなってしまう訳。
また、ちょっとコアな話ですが、誤差が累積していくデータ(たとえばランダムウォークしているデータ)に対して、回帰分析は無力だと 言うことが知られています。
そこで「なんでも回帰分析」ってのはやめ、時系列データの特徴を丁寧にみてやって、時系列モデルを組み立てて行くことになります。
ここではいちいち時系列解析の細かい解説はしません。が、知っとかないとプログラムもできないし予測もできないというくらい大事な 「用語」にかんしては、すこしだけ書いておきます。御参考までに。
自己相
関係数は、過去の値とどれくらい似ているか(あるいは似ていないか)を表したものです。
たとえば、一日前と大きな正の自己相関があれば、「一日前に多ければ、今日も多い」ということになり、2日前と負の自己相関があれば「二日前に多ければ、
今日は少ない」ということになります。
R では acf( データ ) で図示してくれます。
偏自己相関係数は、注目している時以外の要因を無視して計算された自己相関係数です。一日前と今日の関係が知りたくて、一昨日の要因 何か見たくないよってときにはこちらを使います。
Rでは pacf( データ ) で図示してくれます
自己回
帰とは、その名の通り自分のデータと回帰する手法。ようするに、一日前の値を横軸に、一日後のデータを縦軸にしてプロットし、線を引っ張ると言うイメージ
です。
一番簡単な時系列モデル。
Rでは ar( データ , aic=TRUE
) でモデリングしてくれます。 aic=TRUE
がデフォルトなので、AICで勝手にモデル選択(モデル選択を知らない方はこちらをどう
ぞ)をしてくれます。便利。
ただし、step()関数と同じく、全部のモデルを総当たりにAIC計算をしてくれているわけではどうやらなさそう(?)。ようするに、計算量をちょっと
ケチってるんですね。なので、モデル選択してくれた後でも、一応手作業で気になるモデルのAICを計算した方がいいかもしれません。
1期前のデータだけから回帰したARモデルを AR(1) と表します。同様に p期前までの全てのデータを使って今日のデータを予 測するモデルはAR(p) モデルと書きます。
こちらは、ちょっと込み入ったモデルです。というのも、「移動平均」という名前が悪い。これって歴史的な経緯(?)があってこんな名
前になったんですね。
まず、普通の移動平均の説明をします。
4 2 6 10 14 18
こんなデータがあったとします。 これを3区間移動平均をしてみます。すると……
(4+2+6)/3
(2+6+10)/3 (6+10+14)/3 (10+14+18)/3
= 4 6 10 14
となりますこれが移動平均。 ようするに3地点の平均をずらしながらやっていったら3点移動平均になるわけです。
ここで、緑の部分に注目。移動平均する際のたんなる計算過程なのですが・・・・・
4番目の移動平均結果 14
= 10+14+18
÷3
3番目の移動平均結果 10
= 6+10+14
÷3
2番目の移動平均結果
6 = 2+6+10
÷3
1番目の移動平均結果
4 =
4+2+6 ÷3
あおい部分を見てください。4番目の結果に使用されているデータが3番目2番目の移動平均結果にも用いられていることが分かります。
同じ数字が使われている⇒似ている⇒相関がある ということになります。
だから、要するに、「移動平均みたいに、一部同じ数字を使って表されたデータは自己相関を持つ」
ということになるわけです。これを利用したのが移動平均モデル。
ちなみに、ここには書きませんが、式の形は一欠片も移動平均の式に似ていません。「一部同じ数字を使って表されたデータは自己相関を持つ」というとこ
ろだけが一緒です。
MAモデルだと、一日前と大きな自己相関があるとかいう現象や周期性うまくモデリングすることができます。
また、さきほど説明したARモデルとは密接な関係があるのですが、ここでは紹介しません。
3期移動平均モデル(要するに上の例と同じ)はMA(3)と表します。同様にq期移動平均モデルは
MA(q)と表記します。MA(q)はq個同じ数字を使っているわけなので、q期前までとは自己相関があり、q+1期前のデータとは自己相関がありませ
ん。
名前の通り、さっき解説したモデル(ARモデル、MAモデル)を両方突っ込んだモデルです。
データの差分を取ってからARMAを適用したモデルです。
差分を取るのに、なぜ和分? って自分はむかし気になったんですが、要するに視点の問題ですね。もともとの値が累積していって生じた モデルだから和分過程。和分過程を元々のARMA過程に直すためには、差分を取ることになるわけです。
Rでは arima.sim(n=乱数の数, model=list(order=c(p,d,q),ar=c(ARの係数),ma=c(MAの係数)),sd=sqrt(1))
で d 回の和分過程のAR(p)、MA(q) から生じるデータ(以下はARIMA(p、d,q)と記します)を作り出すことができます。シミュレーションに便 利。
実際のモデル推定は
arima(データ , order=c(p,d,q) )
で計算できます orderのあとのp、d、qには、それぞれAR(p)、d 回和分、MA(q)モデルの寄せ集めだということの指定をします。この時のモデルをARIMA(p、d、q)と記します。
ARIMAモデルに、さらに周期的な変動(季節変動とか)を取り入れさせたモデルです。
Rでは arima( データ , order=c(p,d,q),seasonal=list(order=c(pp,dd,qq) )
で 推定可能です。この場合、パラメータの数は p+d+q+pp+dd+qq 個分推定することになります。この数字の組み合わせのパターンだけでも、嫌に なるくらい多くなりますね。AR(1)にするか、AR(2)にするかそれとも2回差分をとるか1回でいいのかも考えないと……って延々数値を突っ込みモデ ルを作り続けるのは大変に面倒です。
そこで、以下の章では便利なパッケージ その名も forecast を紹介します。
forecast パッケージはCRANに登録されているので、
メニューバーのパッケージ⇒パッケージのインストール⇒CRANmirrorから適当にミラーサイトを選ぶ(例えばJAPAN (Tsukuba)など)⇒forecastを探し出してダブルクリック
でインストールできます。
簡単なモデル作成関数として
auto.arima(データ , ic="aic" , trace=T , stepwise=T)
があります名前(auto.arima)の通り、全自動でモデル選択をしてくれます。季節成分があるかも勝手に判断してくれるので、無 茶苦茶便利です。
ic="aic" は、AICを情報量基準として使用してモデル選択をしてねという指定。AIC以外にも、AICc(aiccと入力)とか、BIC(bicと入力)なんかも 使えます。
trace=T はデフォルトではFなのですが、なんとなく有った方が良いような気がします。これをTにしておくと、「こう言うモデ
ルを作ってみてAICを計算してみましたよ」という計算履歴の一覧をずらずらと表示させることができます。
なぜこれが大事かというと、AICの計算もれがあるからですね。履歴を見て気になるところは手作業でAICを計算しましょう。
stepwice=Tは、計算をケチって素早く終わらせなさいよという指示です(デフォルトはT)。これをFにすると、やたら時間が かかってしまいますが、その代り丁寧に計算をしてくれます。私は結構Fにして計算させることが多いですが、どっちの方がいいのかは不明。
やっぱり、自動化されているとはいえ、人間が多少は汗を流した方がよいのかもしれません。
パッ
ケージforecastは、自動SARIMA推定だけでなく、色々な便利な関数が付いています。その中でも、forecast関数(forecastパッ
ケージの中にこんな名前の関数があるのです。ややこしいですね)は大変便利です。将来予測+任意の予測区間の推定が一発で出来上がります。
具体的には
plot(forecast(model,
level = c(,95), h = 50))
で、推定されたモデルの予測値と95%予測区間を50期先まで求めてプロットまでしてくれます。
数字そのものが知りたいならばplotを抜いて。
forecast(model, level = c(,95), h = 50)
ってすれば数字がずらずらと出て来ます。
この関数はlevel=c(95)
の95を変えれば色々な予測幅を計算することができます。たとえば50%予測区間とか。また、 level=c(50,95) とす
れば、50%予測区間も95%予測区間も両方同時に求めることができます。
とても便利なパッケージなので、使ってみては?
参考文献
z。 & nbsp; June−1992
沖本竜義:計量時系列分析、朝倉書店、東京、2010
田中孝文:Rによる時系列分析入門、CAP出版株式会社、東京、2008
RjpWiki アーカイブス 時系列オブジェクト一覧
http://www.okada.jp.org/RWiki/?RjpWiki%20%A5%A2%A1%BC%A5%AB%A5%A4%A5%D6%A5%B9
Package ‘forecast’ 説明書(PDF)
http://cran.r-project.org/web/packages/forecast/forecast.pdf
前のページ ⇒ 平滑化スプラインと加法モデル へ
次のページ ⇒ 時
系列解析 実践編 へ