(-"-)

昨日、一昨日と、全く同じところで寝落ちしました・・・叫び

いつも20時半ごろからとてつもない睡魔が襲ってきて
そこを無理やりお風呂に入ることにより覚醒させ
22時頃から勉強再開～のパターンだったのですが

手を動かさないとだめだということがはっきりしましたあせる

著作権を侵害しない程度に要約したいと思います。

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

今日は残業もあってこの時間スタート。

昨日1章の途中で寝落ちしたので（汗）今日は再度1章から。
1日40ページで10日、数式の理解は諦めてるので最低このくらいのペースで進めたいです。

今日は序章と「経済時系列分析のために」「経済時系列データの構成要素と分解」あたりです。

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

Rでの分析をする前に、一通り概念だけでも把握しておこうと
時系列分析の教科書を読んでいます。

「経済時系列分析のための統計的基礎」として
「確率と数理統計の基礎を持つことが望ましい」と記されているのを読むと
心折れそうになりますが・・・（汗）

ざっと見た限りでも、インターネットであちこち調べた時に見たキーワードが
随所に散らばっているので、理解しやすくなるのではと期待しています。

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

SPSSが作った最適モデルと予測値です。




「ARIMAモデルパラメータ」の解釈の仕方がよくわからず・・・

一旦おいておき、まずは教科書でARとMAからおさらいしたいと思います。

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

単位根をもつデータで回帰分析を行うと見せかけの回帰が発生することがあるとのことで
前期レポート（可処分所得と牛肉・豚肉消費量の回帰分析）はどうだったんだろう？と。


> zai=read.csv("上級財と下級財.csv")
> zai
year shotoku ushi buta
1 2000 474411 26152 21728
2 2001 466003 21157 22396
3 2002 453716 19982 23434
4 2003 440667 21374 21858
5 2004 446288 20918 23362
6 2005 441156 21324 23191
7 2006 441448 20705 23249
8 2007 442504 20868 23923
9 2008 442749 20885 25555
10 2009 427912 20167 24791
11 2010 429967 18964 23957

> adf.test(zai$ushi)

Augmented Dickey-Fuller Test

data: zai$ushi
Dickey-Fuller = 0.49819, Lag order = 2, p-value = 0.99
alternative hypothesis: stationary

Warning message:
In adf.test(zai$ushi) : p-value greater than printed p-value

> adf.test(zai$buta)

Augmented Dickey-Fuller Test

data: zai$buta
Dickey-Fuller = -2.1133, Lag order = 2, p-value = 0.5292
alternative hypothesis: stationary

> adf.test(zai$shotoku)

Augmented Dickey-Fuller Test

data: zai$shotoku
Dickey-Fuller = -2.9715, Lag order = 2, p-value = 0.2023
alternative hypothesis: stationary

> attach(zai)
> tzai=ts(data.frame(shotoku,ushi),start=c(2000),frequency=1)
> tzai.vecm<-ca.jo(tzai,ecdet="none",type="eigen",K=2,spec="longrun")
> summary(tzai.vecm)

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

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

Eigenvalues (lambda):
[1] 0.978205232 0.006646841

Values of teststatistic and critical values of test:

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

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

shotoku.l2 ushi.l2
shotoku.l2 1.00000 1.00000
ushi.l2 -14.70265 14.88364

Weights W:
(This is the loading matrix)

shotoku.l2 ushi.l2
shotoku.d -0.50176474 -0.034212011
ushi.d 0.09580659 -0.002062248


こちらも単位根ありですが、可処分所得と牛肉消費量では共和分ランク１となりました。


> tzai.vec2var<-vec2var(tzai.vecm,r=1)
> print(tzai.vec2var)

Coefficient matrix of lagged endogenous variables:

A1:
shotoku.l1 ushi.l1
shotoku 0.2629252 7.4695378
ushi 0.1215143 0.1326964


A2:
shotoku.l2 ushi.l2
shotoku 0.23531008 -0.09226474
ushi -0.02570774 -0.54130746


Coefficient matrix of deterministic regressor(s).

constant
shotoku 64389.54
ushi -13030.31

> tzai.prd<-predict(tzai.vec2var,n.ahead=5,ci=0.95)
> plot(tzai.prd)


> tzai.irf<-irf(tzai.vec2var,n.ahead=5,ci=0.95)
> plot(tzai.irf)
Hit to see next plot:




一通り結果は出たけど、「可処分所得↑⇒牛肉消費量↑」なのかはわからず・・・
Granger因果性検定をかけようとしたらエラーになりました。変数の型が違うのかも

続きはもう明日にします。

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

階差データはオリジナルの先頭もしくは最終のデータがわからなければ復元できないとは
思いますが、diffinv関数は階差データを復元できるのかと期待したものの、
値を足し合わせていく（階差の逆）関数であることがわかりがっかりしています。

階差データを投入してモデルを作るのは予測のためではなく、
変数間の関係性やショックの影響を調べるためなのかも？

そこへいくと階差をとらずに分析できるというのは大きな魅力ですね。
共和分があってほしかったです。


あと、ACFとPACFを求める関数もありました。

●ACF（インフレ率）
> acf(tmydata[,1])


> acf(tmydata[,1],plot=FALSE)

Autocorrelations of series ‘tmydata[, 1]’, by lag

0 1 2 3 4 5 6 7 8 9 10
1.000 0.727 0.554 0.319 0.234 0.184 0.154 0.132 0.133 0.175 0.197
11 12 13 14 15
0.188 0.087 0.033 -0.064 -0.106

●Ljung-Box検定（インフレ率）
> Box.test(tmydata[,1], lag=15, type="Ljung-Box")

Box-Ljung test

data: tmydata[, 1]
X-squared = 49.653, df = 15, p-value = 1.372e-05

ラグ1～15において「無相関」は棄却。

●ACF（失業率）
> acf(tmydata[,2])


> acf(tmydata[,2],plot=FALSE)

Autocorrelations of series ‘tmydata[, 2]’, by lag

0 1 2 3 4 5 6 7 8 9 10
1.000 0.928 0.811 0.687 0.571 0.483 0.421 0.367 0.293 0.216 0.139
11 12 13 14 15
0.054 -0.033 -0.111 -0.179 -0.235

●Ljung-Box検定（失業率）
> Box.test(tmydata[,2], lag=15, type="Ljung-Box")

Box-Ljung test

data: tmydata[, 2]
X-squared = 129.47, df = 15, p-value < 2.2e-16

ラグ1～15において「無相関」は棄却。

●PACF（インフレ率）
> pacf(tmydata[,1])


> pacf(tmydata[,1],plot=FALSE)

Partial autocorrelations of series ‘tmydata[, 1]’, by lag

1 2 3 4 5 6 7 8 9 10 11
0.727 0.056 -0.217 0.125 0.080 -0.044 0.016 0.076 0.111 0.009 -0.035
12 13 14 15
-0.150 0.031 -0.101 -0.075

●PACF（失業率）
> pacf(tmydata[,2])


> pacf(tmydata[,2],plot=FALSE)

Partial autocorrelations of series ‘tmydata[, 2]’, by lag

1 2 3 4 5 6 7 8 9 10 11
0.928 -0.368 -0.010 -0.011 0.118 0.026 -0.071 -0.220 0.054 -0.031 -0.139
12 13 14 15
-0.111 -0.017 -0.028 -0.004


なんだかちょびっとずつしか進みませんが、地道にやるしかないですよね・・・

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

階差をとって定常になることがわかれば、
gretlで階差をとったデータを用いてモデルの推定をしていたのでできるのかなと思い
とりあえずVARモデルもやってみることにしました。

＞＞＞

●４回の階差をとり、dtmydataへ
> dtmydata<-diff(tmydata,differences=4)

●最適なラグ時数を求める
> VARselect(dtmydata,lag.max=5,type="const")
$selection
AIC(n) HQ(n) SC(n) FPE(n)
5 3 3 3

$criteria
1 2 3 4 5
AIC(n) -7.8877348848 -8.7038693400 -9.0124543099 -8.9609972694 -9.0366540494
HQ(n) -7.8041303583 -8.5645284624 -8.8173770813 -8.7101836899 -8.7301041188
SC(n) -7.5974049145 -8.2199860561 -8.3350177125 -8.0900073585 -7.9721108249
FPE(n) 0.0003760945 0.0001675597 0.0001252448 0.0001362058 0.0001332576

●VARモデルを推定
> dtmydata.var<-VAR(dtmydata,p=5)
> summary(dtmydata.var)

VAR Estimation Results:
=========================
Endogenous variables: inflation, shitsugyo
Deterministic variables: const
Sample size: 26
Log Likelihood: 65.692
Roots of the characteristic polynomial:
0.9592 0.9592 0.8619 0.7874 0.7874 0.7809 0.7809 0.7349 0.7349 0.4464
Call:
VAR(y = dtmydata, p = 5)

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

Estimate Std. Error t value Pr(>|t|)
inflation.l1 -1.927315 0.339514 -5.677 4.39e-05 ***
shitsugyo.l1 -0.009345 0.013917 -0.671 0.51214
inflation.l2 -1.923576 0.563393 -3.414 0.00384 **
shitsugyo.l2 -0.006812 0.022253 -0.306 0.76373
inflation.l3 -1.448192 0.529341 -2.736 0.01532 *
shitsugyo.l3 -0.001956 0.027544 -0.071 0.94432
inflation.l4 -0.948615 0.428317 -2.215 0.04268 *
shitsugyo.l4 -0.008059 0.023203 -0.347 0.73316
inflation.l5 -0.319913 0.209567 -1.527 0.14768
shitsugyo.l5 0.003144 0.014380 0.219 0.82989
const 0.000278 0.004420 0.063 0.95068
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0222 on 15 degrees of freedom
Multiple R-Squared: 0.8736, Adjusted R-squared: 0.7893
F-statistic: 10.37 on 10 and 15 DF, p-value: 4.649e-05

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

Estimate Std. Error t value Pr(>|t|)
inflation.l1 15.63073 8.67003 1.803 0.0915 .
shitsugyo.l1 -1.01870 0.35541 -2.866 0.0118 *
inflation.l2 22.18800 14.38712 1.542 0.1439
shitsugyo.l2 -0.99285 0.56826 -1.747 0.1010
inflation.l3 15.98320 13.51754 1.182 0.2554
shitsugyo.l3 -0.46093 0.70339 -0.655 0.5222
inflation.l4 8.44936 10.93775 0.772 0.4518
shitsugyo.l4 0.14120 0.59252 0.238 0.8149
inflation.l5 3.31836 5.35161 0.620 0.5445
shitsugyo.l5 0.19340 0.36721 0.527 0.6061
const -0.04179 0.11288 -0.370 0.7164
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5669 on 15 degrees of freedom
Multiple R-Squared: 0.8218, Adjusted R-squared: 0.703
F-statistic: 6.918 on 10 and 15 DF, p-value: 0.0004948

Covariance matrix of residuals:
inflation shitsugyo
inflation 0.0004929 -0.009624
shitsugyo -0.0096239 0.321417

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

●インパルス応答関数
> dtmydata.irf<-irf(dtmydata.var,n.ahead=14,ci=0.95)
> plot(dtmydata.irf)
Hit to see next plot:



●グレンジャー因果性検定
> causality(dtmydata.var,cause="inflation")
$Granger

Granger causality H0: inflation do not Granger-cause shitsugyo

data: VAR object dtmydata.var
F-Test = 0.6887, df1 = 5, df2 = 30, p-value = 0.6358

$Instant

H0: No instantaneous causality between: inflation and shitsugyo

data: VAR object dtmydata.var
Chi-squared = 9.5925, df = 1, p-value = 0.001954

●一期先を予測
> dtmydata.pr<-predict(dtmydata.var,n.ahead=1,ci=0.95,dumvar=NULL)
> dtmydata.pr
$inflation
fcst lower upper CI
inflation.fcst -0.003584949 -0.04709803 0.03992813 0.04351308

$shitsugyo
fcst lower upper CI
shitsugyo.fcst -1.762654 -2.873828 -0.651479 1.111175
>
> plot(dtmydata.pr)


＜＜＜

インパルス応答関数とグレンジャー因果性検定はこのままでもいいのかもしれないけど
階差データを使っているので予測値はどうなるんだろう？

VARモデルに投入した配列に予測値の行を追加してdiffinvかなぁと予想するものの
いろいろ調べても方法がわからなくて困っている最中です。

というところで、一旦切りたいと思います。

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

「単位根過程の定義：原系列が非定常過程で、差分系列が定常過程であるもの」
と書かれているのを読んだことで、
非定常過程と単位根過程が頭の中でごちゃまぜになっていたことに気づきました。

改めて教科書を読んでみると、
非定常過程は「単位根過程」「発散過程」に分けられるとありました。

また、diff関数のusageも読んで、こちらも使い方が間違っていたようです。
何回の階差をとったら定常になるかを試してみたのですが・・・

＞＞＞

●１階（インフレ率）
> adf.test(diff(tmydata[,1],differences=1))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 1], differences = 1)
Dickey-Fuller = -3.3976, Lag order = 3, p-value = 0.07454
alternative hypothesis: stationary

●１階（失業率）
> adf.test(diff(tmydata[,2],differences=1))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 2], differences = 1)
Dickey-Fuller = -2.5691, Lag order = 3, p-value = 0.3525
alternative hypothesis: stationary

●２階（インフレ率）
> adf.test(diff(tmydata[,1],differences=2))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 1], differences = 2)
Dickey-Fuller = -4.8082, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(tmydata[, 1], differences = 2)) :
p-value smaller than printed p-value

●２階（失業率）
> adf.test(diff(tmydata[,2],differences=2))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 2], differences = 2)
Dickey-Fuller = -2.7819, Lag order = 3, p-value = 0.2706
alternative hypothesis: stationary

●３階（インフレ率）
> adf.test(diff(tmydata[,1],differences=3))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 1], differences = 3)
Dickey-Fuller = -5.9896, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(tmydata[, 1], differences = 3)) :
p-value smaller than printed p-value

●３階（失業率）
> adf.test(diff(tmydata[,2],differences=3))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 2], differences = 3)
Dickey-Fuller = -3.0798, Lag order = 3, p-value = 0.1562
alternative hypothesis: stationary

●４階（インフレ率）
> adf.test(diff(tmydata[,1],differences=4))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 1], differences = 4)
Dickey-Fuller = -6.4161, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(tmydata[, 1], differences = 4)) :
p-value smaller than printed p-value

●４階（失業率）
> adf.test(diff(tmydata[,2],differences=4))

Augmented Dickey-Fuller Test

data: diff(tmydata[, 2], differences = 4)
Dickey-Fuller = -5.5324, Lag order = 3, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(diff(tmydata[, 2], differences = 4)) :
p-value smaller than printed p-value

＜＜＜

とりあえず全ての要素が定常になる階差は存在し、
それぞれの系列は単位根過程だったようです。（共和分検定の結果には影響しませんが）

計量経済学レポート（１０）

教科書で共和分の定義を読むと、
全ての要素が同じ階差をとると定常になることが必要であるようです。

今行っているのでは、ベストな階差が「2階」でしたが

> dtmydata<-diff(tmydata,lag=2)
> 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

という結果だったので満たせていません。

また、昨日行ったJohansenの共和分検定で「ランク０」という結果になりましたが
ランク０の場合は共和分の関係はないと結論付けられるようです。

VARモデルは不可ということで合ってるかな・・・

また、当時のレポートの結果を見てみると、
SPSSの自動モデル作成機能で作られたモデルは「ARIMA(0,1,1)」ということでした。
時系列分析の教科書を読むとともに、ARIMAモデルについてもやってみたいと思います。

教科書はすごく敷居の高い本と当時は思っていましたが、今必要で調べたいことが
一通りわかるようになっていました。
授業を受け、先生に質問できる環境がある時にちゃんと読んでおけばよかったです。

ちなみに、失業率を被説明変数、インフレ率を説明変数とした時の回帰分析の結果が

> summary(lm(shitsugyo~inflation))

Call:
lm(formula = shitsugyo ~ inflation)

Residuals:
Min 1Q Median 3Q Max
-1.1844 -0.5511 0.2117 0.4849 1.5794

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.9844 0.1429 27.883 < 2e-16 ***
inflation -45.5650 6.3209 -7.209 2.9e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7177 on 33 degrees of freedom
Multiple R-squared: 0.6116, Adjusted R-squared: 0.5998
F-statistic: 51.96 on 1 and 33 DF, p-value: 2.895e-08

となりますが、これが「見せかけの回帰」なのでしょうか・・・