計量経済学レポート（９）

共和分検定の方法と結果はわかったのですが、結果の読み方がわからないので（ランクとか）
まずは教科書で共和分の詳しい概念を読むことにします。


●Johansenの共和分検定（固有値）
> tmydata.vecm<-ca.jo(tmydata,ecdet="none",type="eigen",K=2,spec="longrun")
> summary(tmydata.vecm)

######################
# Johansen-Procedure #
######################

Test type: maximal eigenvalue statistic (lambda max) , with linear trend

Eigenvalues (lambda):
[1] 0.19031127 0.08328484

Values of teststatistic and critical values of test:

test 10pct 5pct 1pct
r <= 1 | 2.87 6.50 8.18 11.65
r = 0 | 6.97 12.91 14.90 19.19

Eigenvectors, normalised to first column:
(These are the cointegration relations)

inflation.l2 shitsugyo.l2
inflation.l2 1.00000000 1.00000000
shitsugyo.l2 0.01024485 0.04666443

Weights W:
(This is the loading matrix)

inflation.l2 shitsugyo.l2
inflation.d -0.4875411 0.0126644
shitsugyo.d 4.5253922 -2.1201067


●Johansenの共和分検定（トレース）
> tmydata.vecm<-ca.jo(tmydata,ecdet="none",type="trace",K=2,spec="longrun")
> summary(tmydata.vecm)

######################
# Johansen-Procedure #
######################

Test type: trace statistic , with linear trend

Eigenvalues (lambda):
[1] 0.19031127 0.08328484

Values of teststatistic and critical values of test:

test 10pct 5pct 1pct
r <= 1 | 2.87 6.50 8.18 11.65
r = 0 | 9.84 15.66 17.95 23.52

Eigenvectors, normalised to first column:
(These are the cointegration relations)

inflation.l2 shitsugyo.l2
inflation.l2 1.00000000 1.00000000
shitsugyo.l2 0.01024485 0.04666443

Weights W:
(This is the loading matrix)

inflation.l2 shitsugyo.l2
inflation.d -0.4875411 0.0126644
shitsugyo.d 4.5253922 -2.1201067

計量経済学レポート（８）

整理のため、最初からやりなおしてみました。

＞＞＞

● ワーキングディレクトリ設定
> getwd()
[1] "C:/Users/・・・・・/Documents"
> setwd("C:\\Users\\・・・・・\\Desktop\\20150509rlec")

● データ読み込み＆時系列データ設定
> mydata=read.csv("完全失業率とインフレ率.csv")
> attach(mydata)
> tmydata=ts(data.frame(inflation,shitsugyo),start=c(1978),frequency=1)

● データ参照
> tmydata
Time Series:
Start = 1978
End = 2012
Frequency = 1
inflation shitsugyo
1978 0.043806647 2.2
1979 0.036179450 2.1
1980 0.078212291 2.0
1981 0.047927461 2.2
1982 0.028430161 2.4
1983 0.018028846 2.6
1984 0.023612751 2.7
1985 0.019607843 2.6
1986 0.006787330 2.8
1987 0.000000000 2.8
1988 0.007865169 2.5
1989 0.022296544 2.3
1990 0.030534351 2.1
1991 0.032804233 2.1
1992 0.017418033 2.2
1993 0.013091641 2.5
1994 0.005964215 2.9
1995 -0.000988142 3.2
1996 0.000989120 3.4
1997 0.018774704 3.4
1998 0.005819593 4.1
1999 -0.002892960 4.7
2000 -0.006769826 4.7
2001 -0.007789679 5.0
2002 -0.008832188 5.4
2003 -0.002970297 5.3
2004 0.000000000 4.7
2005 -0.002979146 4.4
2006 0.002988048 4.1
2007 0.000000000 3.9
2008 0.013902681 4.0
2009 -0.013712047 5.1
2010 -0.006951341 5.1
2011 -0.003000000 4.6
2012 0.000000000 4.3

● varsパッケージインストール＆呼び出し
> install.packages("vars",rep="http://cran.ism.ac.jp")
> library(vars)

● 次数の選択
> VARselect(tmydata,lag.max=5,type="const")
$selection
AIC(n) HQ(n) SC(n) FPE(n)
2 2 1 2

$criteria
1 2 3 4 5
AIC(n) -1.200239e+01 -1.209382e+01 -1.194445e+01 -1.180444e+01 -1.165958e+01
HQ(n) -1.191274e+01 -1.194441e+01 -1.173526e+01 -1.153548e+01 -1.133086e+01
SC(n) -1.172215e+01 -1.162676e+01 -1.129055e+01 -1.096372e+01 -1.063204e+01
FPE(n) 6.137766e-06 5.629185e-06 6.609971e-06 7.761309e-06 9.276880e-06

● VARモデル推定（階差：2）
> tmydata.var<-VAR(tmydata,p=2)
> summary(tmydata.var)

VAR Estimation Results:
=========================
Endogenous variables: inflation, shitsugyo
Deterministic variables: const
Sample size: 33
Log Likelihood: 99.394
Roots of the characteristic polynomial:
0.8043 0.6651 0.4963 0.07686
Call:
VAR(y = tmydata, p = 2)

Estimation results for equation inflation:
==========================================
inflation = inflation.l1 + shitsugyo.l1 + inflation.l2 + shitsugyo.l2 + const

Estimate Std. Error t value Pr(>|t|)
inflation.l1 0.387748 0.216227 1.793 0.0837 .
shitsugyo.l1 -0.014443 0.007599 -1.901 0.0677 .
inflation.l2 0.137376 0.210434 0.653 0.5192
shitsugyo.l2 0.010039 0.007587 1.323 0.1965
const 0.020026 0.013250 1.511 0.1419
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01227 on 28 degrees of freedom
Multiple R-Squared: 0.6237, Adjusted R-squared: 0.5699
F-statistic: 11.6 on 4 and 28 DF, p-value: 1.111e-05

Estimation results for equation shitsugyo:
==========================================
shitsugyo = inflation.l1 + shitsugyo.l1 + inflation.l2 + shitsugyo.l2 + const

Estimate Std. Error t value Pr(>|t|)
inflation.l1 7.9487 5.5850 1.423 0.16572
shitsugyo.l1 1.5011 0.1963 7.648 2.49e-08 ***
inflation.l2 -5.5434 5.4354 -1.020 0.31652
shitsugyo.l2 -0.5536 0.1960 -2.825 0.00862 **
const 0.1884 0.3422 0.551 0.58627
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3168 on 28 degrees of freedom
Multiple R-Squared: 0.9301, Adjusted R-squared: 0.9202
F-statistic: 93.2 on 4 and 28 DF, p-value: 9.248e-16

Covariance matrix of residuals:
inflation shitsugyo
inflation 0.0001504 -0.001891
shitsugyo -0.0018912 0.100373

Correlation matrix of residuals:
inflation shitsugyo
inflation 1.0000 -0.4867
shitsugyo -0.4867 1.0000

● インパルス応答関数（14期先まで・信頼区間95%）
> tmydata.irf<-irf(tmydata.var,n.ahead=14,ci=0.95)
> plot(tmydata.irf)

Hit to see next plot:


● Granger因果性検定
> causality(tmydata.var,cause="inflation")
$Granger

Granger causality H0: inflation do not Granger-cause shitsugyo

data: VAR object tmydata.var
F-Test = 1.0483, df1 = 2, df2 = 56, p-value = 0.3573

$Instant

H0: No instantaneous causality between: inflation and shitsugyo

data: VAR object tmydata.var
Chi-squared = 6.3192, df = 1, p-value = 0.01194

＜＜＜

● tseriesパッケージインストール＆呼び出し
> install.packages("tseries",rep="http://cran.ism.ac.jp")
> library(tseries)

●2回階差を取り、dtmydataに格納
dtmydata<-diff(tmydata,lag=2)

●2階階差をとった時のグラフ
> plot(dtmydata)
Hit to see next plot:


●単位根検定（インフレ率）
> adf.test(dtmydata[,1])

Augmented Dickey-Fuller Test

data: dtmydata[, 1]
Dickey-Fuller = -3.2694, Lag order = 3, p-value = 0.09338
alternative hypothesis: stationary


●単位根検定（失業率）
> adf.test(dtmydata[,2])

Augmented Dickey-Fuller Test

data: dtmydata[, 2]
Dickey-Fuller = -3.4366, Lag order = 3, p-value = 0.06941
alternative hypothesis: stationary

●単位根あり→共和分検定へ

＞＞＞

ここまでがこれまでやったことです。
配列の順番を入れ替えただけで、インパルス応答関数のplotが全然違う形に

計量経済学レポート（７）

VARモデルでインパルス応答関数を導出する際には外生性の高い順にする必要があるそうです。

gretlの計量経済分析の本にも

「gretlでVARモデルを推定したりインパルス応答関数を計算したりする場合には
　経済連鎖の始まりとなりそうな変数から並べるという原則がある」

旨が書かれていたのですが、ツールの利用方法かと思って見落としてました。

思いっきり「失業率→インフレ率」の順でやってました。
いろいろ戻りがありそうです・・・

インパルス応答関数の結果を見ると、
インフレ率が一期先の失業率に影響を与えていそうでした。

ここらへんがモデルの推定にも影響してくるということなのかな

計量経済学レポート（６）

計量経済学の教科書を読んでいるうちに寝落ちしたようです
メガネ踏みつぶされてなくてよかった・・・

ところで昨日の続きですが、
2回の階差がベストのようだけれども、５％水準だと有意になりません。

> dmydata<-diff(mydata,lag=2)
> plot(dmydata)



> adf.test(dmydata[,1])

Augmented Dickey-Fuller Test

data: dmydata[, 1]
Dickey-Fuller = -3.4366, Lag order = 3, p-value = 0.06941
alternative hypothesis: stationary

> adf.test(dmydata[,2])

Augmented Dickey-Fuller Test

data: dmydata[, 2]
Dickey-Fuller = -3.2694, Lag order = 3, p-value = 0.09338
alternative hypothesis: stationary


また、共和分関係があればVECMモデルを推定しVARモデルに変換するという方法が
あるようなのですが、共和分検定の手順がまだわからず調べている最中です。

とっても時間がかかっていますが、中途半端で諦めると気持ち悪さが残りそうなので、
一通りやってみてからR本に移りたいと思います。

一通り終わっても、
「この手順で合っている」という確信は掴めなさそうではありますが・・・

計量経済学レポート（５）

今日は書けることがなさそう・・・

Diffの階差をいくつにしても定常にならない場合には
「あきらめる」という選択肢もあるのでしょうか・・・？

とりあえず、手持ちの本の時系列分析の部分を読んでます。

計量経済学レポート（４）

とりあえずは、定常か非定常かから調べてみる必要があるのかと・・・

> plot(shitsugyo)
> par(new =T)
> plot(inflation)



> adf.test(shitsugyo)

Augmented Dickey-Fuller Test

data: shitsugyo
Dickey-Fuller = -2.156, Lag order = 3, p-value = 0.5125
alternative hypothesis: stationary

> adf.test(inflation)

Augmented Dickey-Fuller Test

data: inflation
Dickey-Fuller = -3.047, Lag order = 3, p-value = 0.1665
alternative hypothesis: stationary


両方単位根をもつようです。

共和分を調べていたらグレンジャー因果性検定というものを見つけたので
やってみました。1回の階差をとっています。


> mydata<-ts(data.frame(shitsugyo,inflation),start=c(1978),frequency=1)
> dmydata<-diff(mydata,lag=1)
> plot(dmydata)



> VARselect(dmydata,lag.max=5)
$selection
AIC(n) HQ(n) SC(n) FPE(n)
1 1 1 1

$criteria
1 2
AIC(n) -1.175734e+01 -1.164552e+01
HQ(n) -1.166874e+01 -1.149786e+01
SC(n) -1.147445e+01 -1.117404e+01
FPE(n) 7.843283e-06 8.819318e-06
3 4
AIC(n) -1.155744e+01 -1.139554e+01
HQ(n) -1.135071e+01 -1.112975e+01
SC(n) -1.089736e+01 -1.054687e+01
FPE(n) 9.752364e-06 1.173178e-05
5
AIC(n) -11.185369389
HQ(n) -10.860512841
SC(n) -10.148110484
FPE(n) 0.000015028

> dmydata.var<-VAR(dmydata,p=1)
> causality(dmydata.var,cause="inflation")
$Granger

Granger causality H0: inflation do not
Granger-cause shitsugyo

data: VAR object dmydata.var
F-Test = 1.6093, df1 = 1, df2 = 60,
p-value = 0.2095


$Instant

H0: No instantaneous causality between:
inflation and shitsugyo

data: VAR object dmydata.var
Chi-squared = 6.5485, df = 1, p-value =
0.0105


インフレは失業率に影響を与えないようです。

逆もやってみました。


> causality(dmydata.var,cause="shitsugyo")
$Granger

Granger causality H0: shitsugyo do not Granger-cause inflation

data: VAR object dmydata.var
F-Test = 2.1685, df1 = 1, df2 = 60, p-value = 0.1461


$Instant

H0: No instantaneous causality between: shitsugyo and inflation

data: VAR object dmydata.var
Chi-squared = 6.5485, df = 1, p-value = 0.0105


失業率もインフレに影響しないとの結果が出ました。（フィリップス曲線は？）

教科書に従って、インパルス応答推定というのもやってみました。


> dmydata.irf<-irf(dmydata.var,n.ahead=14,ci=0.95)
> plot(dmydata.irf)
Hit to see next plot:





この先何をすればよいかわからないので、一旦寝ることにします。

残業

だったので、今日はこれから昨日の続きを開始します。


仕事は楽しく毎日充実しているし、一緒に外食できる友達もいるし
とても幸せなはずなのですが・・・

仕事や勉強から離れる隙間の時間、
なんだかさみしく満たされない気持ちになることが多くなってきました。

そのうち南国から離れることも考えようかなぁ～

なんてことを漠然と考えたりしています。
電カル対応があるのでもちろん今は無理ですが。。。

計量経済学レポート（３）

「VARモデルに含まれる全ての時系列に対して単位根検定を実施する必要がある」

って、どこかのサイトに書かれていました

しかし、「共和分というものが認められれば単位根があってもよい」（？）
というようなことが教科書には書かれています。

ところが
「単位根過程を含むVARモデルは仮に共和分関係を含んでいる場合、
共和分関係を修正したモデルを新たに推定してVARモデルに変換する必要がある」
と書かれているのも読んで、もう何をどうしたらよいかわからなくなってきました。

今日はネットばかりしていたので、明日は本を読む日にしたいと思います。

本は来たけど

復習後に開始したいと思います。

10月の情報セキュリティスペシャリスト試験リベンジのため、Rの勉強は６月いっぱいまでとし、情報処理試験が終わるまで一旦休止する予定です。

時間が限られているので、６月はプライベートの時間がRで埋まるくらいに勉強したいと思ってます。

計量経済学レポート（２）

今度は後期のレポートをRでやってみたいと思います。

任意の時系列データを分析せよというもので（前期は「回帰分析」指定だったかも）
当時経済政策の授業でやっていたフィリップス曲線の理論を用いることを思いつき、
1978～2012年の分析結果から、2013年のインフレ率を「2%」とした時の完全失業率を
予測するというものをやって提出していました。

ただ、当時のレポートを読むとふざけたことに
「SPSSの時系列分析ツールに含まれる自動モデル作成機能を利用し」
とか書かれています

RではVARモデルでやってみました。

＞＞＞

（１）パッケージのインストール＆呼び出し

> install.packages("vars",rep="http://cran.ism.ac.jp")
> library(vars)


（２）データの読み込み

> mydata=read.csv("完全失業率とインフレ率.csv")
> mydata
year shitsugyo inflation
1 1978 2.2 0.043806647
2 1979 2.1 0.036179450
3 1980 2.0 0.078212291
　　　　　　　・
　　　　　　　・
　　　　　　　・


（３）時系列オブジェクトへの変換

> attach(mydata)
> tsdata<-ts(data.frame(shitsugyo,inflation),start=c(1978),frequency=1)
> tsdata
Time Series:
Start = 1978
End = 2012
Frequency = 1
shitsugyo inflation
1978 2.2 0.043806647
1979 2.1 0.036179450
1980 2.0 0.078212291
1981 2.2 0.047927461
1982 2.4 0.028430161
1983 2.6 0.018028846
1984 2.7 0.023612751
1985 2.6 0.019607843
1986 2.8 0.006787330
1987 2.8 0.000000000
1988 2.5 0.007865169
1989 2.3 0.022296544
1990 2.1 0.030534351
1991 2.1 0.032804233
1992 2.2 0.017418033
1993 2.5 0.013091641
1994 2.9 0.005964215
1995 3.2 -0.000988142
1996 3.4 0.000989120
1997 3.4 0.018774704
1998 4.1 0.005819593
1999 4.7 -0.002892960
2000 4.7 -0.006769826
2001 5.0 -0.007789679
2002 5.4 -0.008832188
2003 5.3 -0.002970297
2004 4.7 0.000000000
2005 4.4 -0.002979146
2006 4.1 0.002988048
2007 3.9 0.000000000
2008 4.0 0.013902681
2009 5.1 -0.013712047
2010 5.1 -0.006951341
2011 4.6 -0.003000000
2012 4.3 0.000000000


（４）最適な次数を調べる

> VARselect(tsdata,lag.max=5,type="const")
$selection
AIC(n) HQ(n) SC(n) FPE(n)
2 2 1 2

$criteria
1 2 3 4 5
AIC(n) -1.200239e+01 -1.209382e+01 -1.194445e+01 -1.180444e+01 -1.165958e+01
HQ(n) -1.191274e+01 -1.194441e+01 -1.173526e+01 -1.153548e+01 -1.133086e+01
SC(n) -1.172215e+01 -1.162676e+01 -1.129055e+01 -1.096372e+01 -1.063204e+01
FPE(n) 6.137766e-06 5.629185e-06 6.609971e-06 7.761309e-06 9.276880e-06


（５）VARモデルを推定

> tsdata.var<-VAR(tsdata,p=2,type="const")
> summary(tsdata.var)

VAR Estimation Results:
=========================
Endogenous variables: shitsugyo, inflation
Deterministic variables: const
Sample size: 33
Log Likelihood: 99.394
Roots of the characteristic polynomial:
0.8043 0.6651 0.4963 0.07686
Call:
VAR(y = tsdata, p = 2, type = "const")


Estimation results for equation shitsugyo:
==========================================
shitsugyo = shitsugyo.l1 + inflation.l1 + shitsugyo.l2 + inflation.l2 + const

Estimate Std. Error t value Pr(>|t|)
shitsugyo.l1 1.5011 0.1963 7.648 2.49e-08 ***
inflation.l1 7.9487 5.5850 1.423 0.16572
shitsugyo.l2 -0.5536 0.1960 -2.825 0.00862 **
inflation.l2 -5.5434 5.4354 -1.020 0.31652
const 0.1884 0.3422 0.551 0.58627
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.3168 on 28 degrees of freedom
Multiple R-Squared: 0.9301, Adjusted R-squared: 0.9202
F-statistic: 93.2 on 4 and 28 DF, p-value: 9.248e-16


Estimation results for equation inflation:
==========================================
inflation = shitsugyo.l1 + inflation.l1 + shitsugyo.l2 + inflation.l2 + const

Estimate Std. Error t value Pr(>|t|)
shitsugyo.l1 -0.014443 0.007599 -1.901 0.0677 .
inflation.l1 0.387748 0.216227 1.793 0.0837 .
shitsugyo.l2 0.010039 0.007587 1.323 0.1965
inflation.l2 0.137376 0.210434 0.653 0.5192
const 0.020026 0.013250 1.511 0.1419
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.01227 on 28 degrees of freedom
Multiple R-Squared: 0.6237, Adjusted R-squared: 0.5699
F-statistic: 11.6 on 4 and 28 DF, p-value: 1.111e-05



Covariance matrix of residuals:
shitsugyo inflation
shitsugyo 0.100373 -0.0018912
inflation -0.001891 0.0001504

Correlation matrix of residuals:
shitsugyo inflation
shitsugyo 1.0000 -0.4867
inflation -0.4867 1.0000


（６）一期先を予測（95%信頼区間）

> predict(tsdata.var,n.ahead=1,ci=0.95,dumvar=NULL)
$shitsugyo
fcst lower upper CI
shitsugyo.fcst 4.112916 3.491968 4.733865 0.6209484

$inflation
fcst lower upper CI
inflation.fcst 0.003688914 -0.02035151 0.02772934 0.02404043


（７）予測結果をプロット

> plot(predict(tsdata.var,n.ahead=1,ci=0.95,dumvar=NULL))



＜＜＜

SPSSの結果と大きくずれてはいないものの一致はしていません。
（予測4.03、上4,59、下3.52）

ネット上のサンプルコード丸写しで取りあえずやってみたけど
明日教科書を読んで咀嚼したいと思います。