R、Pythonの勉強の過程をメモっていくノート。2018-05-04RでVARモデル&インパルス応答関数を求めるとき、「ショックの値」をどうやって出すのか?RでVARモデルを推定してインパルス応答関数を出す時に、インパルス応答関数が対応しているところの「ショック」の大きさが幾らなのかを、どこから得たらいいのかという疑問がわきました。いう図を見せられたとして、まぁ増えるのか減るのかさえ分かればいい場合が多いかもしれませんが、タテ軸の数字はレスポンス側の変数の値だから意味が分かるものの、これがインパルス側の変数の当期の値に「どれだけのショック」を与えた場合に生まれる変動なのか数字で言ってくれ、って頼まれたらどう答えるのかなと。Pfaff氏)にメールで「ショックの大きさは、撹乱項の標準偏差の推定値ってことで合ってますか?」と聞いてみました。と打って俺のPsiのコードを見とけ!」と言われました。ここでいうPsiは、https://cran.r-project.org/web/packages/vars/vars.pdfこの仕様書の28ページの上の方に載ってる式のΨのことで、この行列の中身が直交化インパルス応答に対応しています。これはVARのVMA表現ですね。return(Psi)}ここでsigma.uとなっているものが、残差の分散共分散行列ですね!ということは、この行列と同じものを自分で求めて、対角要素の平方根を取ればショックの値になるはずです。これは、VARモデルの推定結果であるx(varestクラスオブジェクト)の中に入っている残差の行列resid(x)を取り出してクロス積をとって自由度で割ったものです。+1は定数項の分になります。以下、念のため実際にサンプルデータで推定してみてやり方を確認しておくことにします。データを読み込みます。変数が4個ありますが、労働生産性と実質賃金だけ使おうと思います。どんなデータか描画しておきます。実質賃金では、VARモデルを推定していきます。簡単にするため、労働生産性(prod)と実質賃金(rw)という2つの変数だけでやってみます。VAR(3)の推定が終わりました。いう方向の、直交インパルス応答関数を求めます。あわせてグラフの描画もしておきます。plot(irf.canada)これでつまり、当期(第1期)のprod(労働生産性)に撹乱項の1標準偏差分のショックが加えられた時の、各期のrw(実質賃金)の増減が得られました。問題は「その労働生産性の撹乱項1標準偏差分ってどんだけやねん」ということですが、これは次のようにして計算できます。ここまでで、ショックの値は計算できました。さっきのインパルス応答関数と合わせていうと、当期の労働生産性の撹乱項に0.684086のショックを加えると、当期の実質賃金には0.003789932、2期目には-0.083499262...の影響があるということが分かったわけです。(実際は先に単位根とか共和分とかを見ないといけませんがここでは分析自体が目的ではないので省略。)ちなみにこの撹乱項の標準偏差の値、どっかに入ってるのではないかと思って探してみたら、いちいち計算しなくても、以下のようにすれば一発で取り出すことができました。実務上はこれを使えばいいと思います*2。error'として載っています*3。念のためですが、さっきのsigma.uの対角要素は分散、こっちのsigmaは標準偏差の推定値です。沖本(2010).経済・ファイナンスデータの計量時系列分析,朝倉書店.に書いてある方法に照らしながら手計算でインパルス応答を出してみて、それが{vars}パッケージのirf()関数の結果と確認するかを見ておきます。後述の通り、使われているショックの値が合ってるかどうかを確認するだけなら、1年後、2年後...の逐次計算をしなくても、当期のresponseだけ見れば分かります。沖本(2010)のp.87-90あたりに書いてある「三角分解」(ほかの呼び方もあるようですが)の行列Aは、以下のように作れますね。matrix(c(1,a,0,1),2,2)で、1つ目の変数から2つ目の変数*4へのインパルス応答のうち、「当期の値」(今回の例でいうと1期目のrwのレスポンスの値)については、「2つ目の変数の撹乱項のうち、1つ目の変数から影響を受けてる部分」がその値に一致すると思います。沖本(2010)のp.89の、この部分に「0.3u_1t」ってのがありますが、要するに2つ目の変数の撹乱項ってのは、2つ目の変数に固有の撹乱項と、1つ目の変数の撹乱項を0.3倍した成分に分かれるということですね。これを分けたのが直行化ってやつです。インパルス応答関数とVARモデルの定義から言って、2つ目の変数の1期目のレスポンスには、1つ目の変数の1期目に与えられたショックを0.3倍(今回の場合でいうとa倍)した値だけが入ることになるので、これだけなら手計算がすぐできます。2期目以降の手計算は面倒ですが。今回の例でいうと、以下のようになります。この通り、さっきの「0.003789932」と一致してるので、合ってると思います。いつも忘れてしまうのですが、非直交化インパルス応答関数であれば、インパルス側の当期の値にショックがあっても、レスポンス側の当期の値は変わりません。直交化インパルス応答関数は、インパルスとレスポンスの攪乱項同士に相関があることを仮定するものなので、当期のインパルス側変数にショックが入ると、当期のレスポンス側変数にもちょっとだけ影響があります。Rで出力されるグラフの、「1」というのは当期のことで、「10」は第10期、つまり9期先ということですね。*1:VMA表現についての解説はハミルトン『時系列解析』の10-11章を見るのがいいと思います。*2:個人的には、撹乱項の標準偏差を推定する上で自由度をいくらにするか自信がなかったので、確認できてよかったですw*3:eとかrwの推定値の標準誤差だから、撹乱項の標準偏差のことになります。*4:沖本(2010)のp.90でいう再帰性、つまり変数の順序性があるので、1つ目、2つ目という意識は意外と重要ですね。最近、パネルデータを扱う...2019-08-08Rで連立方程式を解く練習(例:超小型の日本経済マクロ計量モデル)先日、『gretlで計量経済分析』という本に載っている二段階最小...2017-03-17フォルダの差分同期をExcel(VBA)で行うVBAの勉強エントリです。自分がデータ分析とかをやる場面では...2017-03-07メールデータ解析のため、Outlookの分類フォルダをまたいでメッセージを一括テキスト変換Outlookのマクロ(VBA)に関するエントリです。R、Pythonの勉強の過程をメモっていくノートです。初心者なのであまり信用なさらぬよう。間違いについてはコメント等でご指摘頂ければ幸いです。(1)月別アーカイブつぶやき@statsbeginnerさんのツイートはてなブログをはじめよう!midnightseminarさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。