データの読み込み
capital <- read.csv("data/capitalkuma.csv", header=T)※以下で記載されているアウトライン番号やページ数は論文中のもの
変数の作成
capital$K1_N14 <- capital$K1/capital$N14*100
capital$A1_K1 <- capital$A1/capital$K1*100
capital$S1_C1 <- capital$S1/capital$C1*100
capital$U_L <- capital$U/capital$L*100
capital$I_CP <- (capital$I/capital$I[capital$YEAR==2015])/(capital$CP/capital$CP[capital$YEAR==2015])*100 # 基準値を2015年に
capital$P_N <- capital$P/capital$N #「人口千人につき」という形
capital$M2029_N <- capital$M2029/capital$N*100
capital$E1_N <- capital$E1/(capital$N*1000)*100 # Nが千人なので1000倍しているデータの年代を松村・竹内(1990)にそろえる
matsu <- subset(capital, 1953<=YEAR & YEAR<=1987)library(RcmdrMisc)
matsu_sum <- numSummary(matsu[,c("K1_N14","A1_K1","S1_C1","I_CP","U_L","P_N","M2029_N","E1_N"),drop=F],statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.5,1))
matsu_sum1 <- cbind(matsu_sum$table,matsu_sum$n)
print(round(matsu_sum1,3))## mean sd 0% 50% 100%
## K1_N14 2.593 0.893 1.512 2.257 4.634 35
## A1_K1 97.451 0.671 95.906 97.569 98.710 35
## S1_C1 0.970 0.598 0.221 0.830 2.545 35
## I_CP 65.700 23.217 30.734 67.109 101.128 35
## U_L 1.852 0.556 1.118 1.880 2.844 35
## P_N 14.830 3.313 10.356 12.961 22.084 35
## M2029_N 8.411 0.941 6.607 8.891 9.477 35
## E1_N 1.410 0.516 0.587 1.625 1.963 35matsu_ols <- lm(K1_N14~A1_K1+S1_C1+I_CP+U_L+P_N+M2029_N+E1_N,data=matsu)
library(memisc) #表示を整える
matsu_model <- mtable("1"=matsu_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(matsu_model)##
## Calls:
## 1: lm(formula = K1_N14 ~ A1_K1 + S1_C1 + I_CP + U_L + P_N + M2029_N +
## E1_N, data = matsu)
##
## ==============================
## (Intercept) -7.068*
## (-1.775)
## A1_K1 0.037
## (0.895)
## S1_C1 0.086
## (1.169)
## I_CP -0.000
## (-0.023)
## U_L 0.488***
## (5.944)
## P_N 0.158***
## (4.285)
## M2029_N 0.371***
## (4.863)
## E1_N -0.290
## (-1.005)
## ------------------------------
## adj. R-squared 0.974
## N 35
## ==============================
## Significance:
## *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
library(lmtest)
dwtest(matsu_ols)##
## Durbin-Watson test
##
## data: matsu_ols
## DW = 1.4282, p-value = 0.004377
## alternative hypothesis: true autocorrelation is greater than 0変数の準備（対数変換等）
matsu$S1_C1_LAG <- c(NA,matsu$S1_C1[1:length(matsu$S1_C1)-1])
matsu$LK1_N14 <- log(matsu$K1_N14)
matsu$LA1_K1 <- log(matsu$A1_K1)
matsu$LS1_C1 <- log(matsu$S1_C1)
matsu$LI_CP <- log(matsu$I_CP)
matsu$LU_L <- log(matsu$U_L)
matsu$LP_N <- log(matsu$P_N)
matsu$LM2029_N <- log(matsu$M2029_N)
matsu$LE1_N <- log(matsu$E1_N)
matsu$LS1_C1_LAG <- log(matsu$S1_C1_LAG)表1-5のうち「(2)線形ラグ」
matsu_lag_ols <- lm(K1_N14~A1_K1+S1_C1_LAG+I_CP+U_L+P_N+M2029_N+E1_N,data=matsu)
matsu_lag_model <- mtable("1"=matsu_lag_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(matsu_lag_model)##
## Calls:
## 1: lm(formula = K1_N14 ~ A1_K1 + S1_C1_LAG + I_CP + U_L + P_N +
## M2029_N + E1_N, data = matsu)
##
## ==============================
## (Intercept) -6.046
## (-1.540)
## A1_K1 0.043
## (1.041)
## S1_C1_LAG -0.101
## (-1.271)
## I_CP -0.011
## (-0.999)
## U_L 0.592***
## (6.585)
## P_N 0.116**
## (2.666)
## M2029_N 0.336***
## (4.406)
## E1_N -0.296
## (-1.032)
## ------------------------------
## adj. R-squared 0.972
## N 34
## ==============================
## Significance:
## *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1表1-5のうち「(1)対数」と「(3)対数ラグ」
matsu_log_ols <- lm(LK1_N14~LA1_K1+LS1_C1+LI_CP+LU_L+LP_N+LM2029_N+LE1_N,data=matsu)
matsu_laglog_ols <- lm(LK1_N14~LA1_K1+LS1_C1_LAG+LI_CP+LU_L+LP_N+LM2029_N+LE1_N,data=matsu)
matsu_log_model <- mtable("1"=matsu_log_ols,"2"=matsu_laglog_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(matsu_log_model)##
## Calls:
## 1: lm(formula = LK1_N14 ~ LA1_K1 + LS1_C1 + LI_CP + LU_L + LP_N +
## LM2029_N + LE1_N, data = matsu)
## 2: lm(formula = LK1_N14 ~ LA1_K1 + LS1_C1_LAG + LI_CP + LU_L + LP_N +
## LM2029_N + LE1_N, data = matsu)
##
## ==========================================
## 1 2
## ------------------------------------------
## (Intercept) -2.771 -1.777
## (-0.505) (-0.332)
## LA1_K1 0.065 -0.003
## (0.053) (-0.003)
## LS1_C1 0.000
## (0.007)
## LI_CP -0.107 -0.201
## (-0.541) (-1.008)
## LU_L 0.209*** 0.211***
## (5.275) (5.433)
## LP_N 0.589*** 0.525**
## (2.854) (2.544)
## LM2029_N 1.025*** 0.962***
## (5.941) (5.644)
## LE1_N -0.216* -0.205*
## (-1.943) (-1.881)
## LS1_C1_LAG -0.028
## (-1.436)
## ------------------------------------------
## adj. R-squared 0.982 0.982
## N 35 34
## ==========================================
## Significance: *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1注21：分散不均一性の検定 (p.353)
bptest(matsu_ols)##
## studentized Breusch-Pagan test
##
## data: matsu_ols
## BP = 13.604, df = 7, p-value = 0.0587matsu_cor <- cor(matsu[,c("A1_K1","S1_C1","I_CP","U_L","P_N","M2029_N","E1_N")])
print(matsu_cor,digits=3)## A1_K1 S1_C1 I_CP U_L P_N M2029_N E1_N
## A1_K1 1.0000 0.2739 -0.372 -0.0308 0.430 0.204 -0.401
## S1_C1 0.2739 1.0000 -0.735 0.0549 0.778 0.299 -0.780
## I_CP -0.3724 -0.7348 1.000 0.3871 -0.943 -0.717 0.960
## U_L -0.0308 0.0549 0.387 1.0000 -0.145 -0.749 0.200
## P_N 0.4295 0.7783 -0.943 -0.1448 1.000 0.503 -0.962
## M2029_N 0.2042 0.2989 -0.717 -0.7493 0.503 1.000 -0.535
## E1_N -0.4008 -0.7799 0.960 0.2003 -0.962 -0.535 1.000library(car)
matsu_ols_vif <- vif(matsu_ols)
print(matsu_ols_vif)## A1_K1 S1_C1 I_CP U_L P_N M2029_N E1_N
## 1.295816 3.195856 94.295652 3.476174 24.941617 8.607538 36.825657トレンドと定数項
library(CADFtest)
CADFtest(matsu$K1_N14, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$K1_N14
## ADF(2) = -1.4878, p-value = 0.8114
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1526712CADFtest(matsu$A1_K1, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$A1_K1
## ADF(0) = -4.9915, p-value = 0.001881
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.9294526CADFtest(matsu$S1_C1, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$S1_C1
## ADF(3) = -1.8899, p-value = 0.6347
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.446498CADFtest(matsu$I_CP, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$I_CP
## ADF(0) = -1.1427, p-value = 0.9043
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.114328CADFtest(matsu$U_L, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$U_L
## ADF(2) = -1.668, p-value = 0.7404
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1239CADFtest(matsu$P_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$P_N
## ADF(2) = -1.1861, p-value = 0.8954
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.08833341CADFtest(matsu$M2029_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$M2029_N
## ADF(2) = -1.5467, p-value = 0.7899
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.07461354CADFtest(matsu$E1_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$E1_N
## ADF(2) = -1.9471, p-value = 0.6054
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0596375定数項
CADFtest(matsu$K1_N14, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$K1_N14
## ADF(0) = -1.3596, p-value = 0.5883
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.04780755CADFtest(matsu$A1_K1, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$A1_K1
## ADF(0) = -4.6919, p-value = 0.0007476
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.8616541CADFtest(matsu$S1_C1, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$S1_C1
## ADF(2) = -3.1176, p-value = 0.0359
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.4421995CADFtest(matsu$I_CP, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$I_CP
## ADF(0) = -0.88366, p-value = 0.7795
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01200555CADFtest(matsu$U_L, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$U_L
## ADF(3) = -0.73543, p-value = 0.8226
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.05199654CADFtest(matsu$P_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$P_N
## ADF(3) = -0.37319, p-value = 0.9016
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01129573CADFtest(matsu$M2029_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$M2029_N
## ADF(2) = -0.54812, p-value = 0.8677
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01864944CADFtest(matsu$E1_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$E1_N
## ADF(2) = -1.6329, p-value = 0.4539
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01606601どちらもなし
CADFtest(matsu$K1_N14, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$K1_N14
## ADF(0) = -3.3297, p-value = 0.00165
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.02984538CADFtest(matsu$A1_K1, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$A1_K1
## ADF(3) = -0.39817, p-value = 0.5318
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0005409986CADFtest(matsu$S1_C1, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$S1_C1
## ADF(1) = -2.4983, p-value = 0.01434
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1749909CADFtest(matsu$I_CP, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$I_CP
## ADF(3) = 1.2018, p-value = 0.9375
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.0108353CADFtest(matsu$U_L, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$U_L
## ADF(3) = 0.49158, p-value = 0.8156
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.008961667CADFtest(matsu$P_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$P_N
## ADF(3) = -1.7247, p-value = 0.08
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01112568CADFtest(matsu$M2029_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$M2029_N
## ADF(1) = -0.97612, p-value = 0.287
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0033662CADFtest(matsu$E1_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: matsu$E1_N
## ADF(1) = 0.65028, p-value = 0.8511
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.002364882トレンドと定数項
library(urca)
summary(ur.kpss(matsu$K1_N14,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.2335
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$A1_K1,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.0765
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$S1_C1,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.2057
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$I_CP,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.1182
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$U_L,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.2334
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$P_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.1961
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$M2029_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.2161
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(matsu$E1_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.1928
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216定数項
summary(ur.kpss(matsu$K1_N14,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.9135
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$A1_K1,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.4367
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$S1_C1,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.7456
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$I_CP,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.9717
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$U_L,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.37
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$P_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.9118
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$M2029_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.6436
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(matsu$E1_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 3 lags.
##
## Value of test-statistic is: 0.9
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739変数の作成
capital$K1_N <- capital$K1/capital$N*100
capital$K2_N <- capital$K2/capital$N*100
capital$K3_N <- capital$K3/capital$N*100
capital$A2_K2 <- capital$A2/capital$K2*100
capital$A3_K3 <- capital$A3/capital$K3*100
capital$C1_A1 <- capital$C1/capital$A1*100
capital$C2_A2 <- capital$C2/capital$A2*100
capital$C3_A3 <- capital$C3/capital$A3*100
capital$S1_C1 <- capital$S1/capital$C1*100
capital$S2_C2 <- capital$S2/capital$C2*100
capital$S3_C3 <- capital$S3/capital$C3*100
capital$X_C1 <- capital$X/capital$C1*100
capital$X_C2 <- capital$X/capital$C2*100
capital$X_C3 <- capital$X/capital$C3*100
capital$G1_N <- capital$G1/capital$N*100
capital$G2_N <- capital$G2/capital$N*100
capital$F1_N <- capital$F1/(capital$N*1000)*100
capital$F2_N <- capital$F2/(capital$N*1000)*100
capital$F3_N <- capital$F3/(capital$N/10)*100
capital$M1524_N <- capital$M1524/capital$N*100
capital$E2_N <- capital$E2/(capital$N*1000)*100
capital$Y1972 <- ifelse(capital$YEAR<1972,0,1) #沖縄返還データの年代を秋葉(1993)にそろえる
akiba <- subset(capital, 1960<=YEAR & YEAR<=1986)その他の変数の作成
akiba$TREND=1:nrow(akiba) # トレンド変数library(RcmdrMisc)
akiba_sum <- numSummary(akiba[,c("K1_N","K2_N","K3_N","A1_K1","A2_K2","A3_K3","C1_A1","C2_A2","C3_A3","S1_C1","S2_C2","S3_C3","X_C1","X_C2","X_C3","U_L","W","G1_N","G2_N","F1_N","F2_N","F3_N","M1524_N","E2_N"),drop=F],statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.5,1))
akiba_sum1 <- cbind(akiba_sum$table,akiba_sum$n)
print(round(akiba_sum1,3))## mean sd 0% 50% 100%
## K1_N 1.734 0.355 1.286 1.669 2.510 27
## K2_N 1.925 0.416 1.377 1.874 2.835 27
## K3_N 1.666 0.328 1.247 1.630 2.365 27
## A1_K1 97.324 0.640 95.906 97.295 98.596 27
## A2_K2 96.846 0.615 95.405 96.939 97.938 27
## A3_K3 97.702 0.604 96.458 97.686 98.930 27
## C1_A1 62.721 3.926 54.425 63.310 73.104 27
## C2_A2 55.543 3.684 47.233 55.660 64.401 27
## C3_A3 61.340 4.209 51.829 61.889 72.660 27
## S1_C1 0.751 0.389 0.221 0.763 1.988 27
## S2_C2 0.284 0.157 0.000 0.296 0.561 27
## S3_C3 0.293 0.162 0.000 0.306 0.579 27
## X_C1 0.738 0.817 0.000 0.319 2.648 27
## X_C2 0.761 0.849 0.000 0.333 2.778 27
## X_C3 0.785 0.876 0.000 0.343 2.893 27
## U_L 1.736 0.562 1.118 1.447 2.774 27
## W 57.911 18.979 29.000 68.200 80.400 27
## G1_N 190.131 66.572 76.733 206.660 290.026 27
## G2_N 122.023 86.703 17.138 103.485 275.706 27
## F1_N 8.819 3.284 5.156 7.885 15.566 27
## F2_N 8.717 3.151 5.156 7.885 14.261 27
## F3_N 7.374 3.210 3.698 6.053 13.627 27
## M1524_N 8.492 1.284 6.888 8.541 10.265 27
## E2_N 1.568 0.395 0.760 1.754 1.947 27akiba1_S1_ols <- lm(K1_N~A1_K1+C1_A1+S1_C1+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba2_S2_ols <- lm(K2_N~A2_K2+C2_A2+S2_C2+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba3_S3_ols <- lm(K3_N~A3_K3+C3_A3+S3_C3+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
library(memisc)
akiba_S_model <- mtable("1"=akiba1_S1_ols,"2"=akiba2_S2_ols,"3"=akiba3_S3_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_S_model)##
## Calls:
## 1: lm(formula = K1_N ~ A1_K1 + C1_A1 + S1_C1 + U_L + W + G1_N +
## F1_N + M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = K2_N ~ A2_K2 + C2_A2 + S2_C2 + U_L + W + G1_N +
## F1_N + M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 3: lm(formula = K3_N ~ A3_K3 + C3_A3 + S3_C3 + U_L + W + G1_N +
## F1_N + M1524_N + E2_N + Y1972 + TREND, data = akiba)
##
## ===================================================
## 1 2 3
## ---------- ---------- ----------
## K1_N K2_N K3_N
## ---------------------------------------------------
## (Intercept) 1.313 1.060 1.356
## (0.431) (0.257) (0.356)
## A1_K1 -0.007
## (-0.296)
## C1_A1 -0.007*
## (-1.981)
## S1_C1 0.068
## (1.540)
## U_L 0.370** 0.359** 0.327**
## (2.684) (2.534) (2.415)
## W 0.013 0.015 0.013
## (1.053) (0.983) (1.005)
## G1_N -0.003 -0.002 -0.002
## (-0.561) (-0.386) (-0.361)
## F1_N 0.025 0.063 0.028
## (0.284) (0.691) (0.306)
## M1524_N 0.133** 0.091 0.100*
## (2.379) (1.486) (1.919)
## E2_N -0.079 -0.027 -0.111
## (-0.216) (-0.068) (-0.285)
## Y1972 0.123 0.065 0.065
## (0.999) (0.461) (0.488)
## TREND -0.046 -0.052 -0.049
## (-0.908) (-0.959) (-0.951)
## A2_K2 -0.002
## (-0.073)
## C2_A2 -0.011**
## (-2.561)
## S2_C2 0.139
## (1.394)
## A3_K3 -0.004
## (-0.156)
## C3_A3 -0.009**
## (-2.420)
## S3_C3 0.109
## (1.190)
## ---------------------------------------------------
## adj. R-squared 0.970 0.977 0.965
## N 27 27 27
## ===================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービンワトソン比
library(lmtest)
dwtest(akiba1_S1_ols)##
## Durbin-Watson test
##
## data: akiba1_S1_ols
## DW = 2.0519, p-value = 0.101
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba2_S2_ols)##
## Durbin-Watson test
##
## data: akiba2_S2_ols
## DW = 2.1449, p-value = 0.1566
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba3_S3_ols)##
## Durbin-Watson test
##
## data: akiba3_S3_ols
## DW = 1.9525, p-value = 0.05968
## alternative hypothesis: true autocorrelation is greater than 0akiba1_X_ols <- lm(K1_N~A1_K1+C1_A1+X_C1+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba2_X_ols <- lm(K2_N~A2_K2+C2_A2+X_C2+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba3_X_ols <- lm(K3_N~A3_K3+C3_A3+X_C3+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba_X_model <- mtable("1"=akiba1_X_ols,"2"=akiba2_X_ols,"3"=akiba3_X_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_X_model)##
## Calls:
## 1: lm(formula = K1_N ~ A1_K1 + C1_A1 + X_C1 + U_L + W + G1_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + U_L + W + G1_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 3: lm(formula = K3_N ~ A3_K3 + C3_A3 + X_C3 + U_L + W + G1_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
##
## ===================================================
## 1 2 3
## ---------- ---------- ----------
## K1_N K2_N K3_N
## ---------------------------------------------------
## (Intercept) 2.356 1.341 3.063
## (0.779) (0.339) (0.805)
## A1_K1 -0.023
## (-0.909)
## C1_A1 -0.009**
## (-2.563)
## X_C1 -0.044*
## (-1.852)
## U_L 0.331** 0.344** 0.321**
## (2.603) (2.632) (2.545)
## W 0.012 0.011 0.010
## (0.975) (0.772) (0.817)
## G1_N -0.003 -0.002 -0.001
## (-0.518) (-0.302) (-0.159)
## F1_N 0.094 0.114 0.055
## (1.119) (1.292) (0.628)
## M1524_N 0.091* 0.068 0.080
## (1.794) (1.145) (1.617)
## E2_N 0.165 0.231 0.027
## (0.460) (0.614) (0.074)
## Y1972 0.089 0.090 0.080
## (0.726) (0.692) (0.659)
## TREND -0.034 -0.045 -0.052
## (-0.720) (-0.891) (-1.048)
## A2_K2 -0.011
## (-0.363)
## C2_A2 -0.010**
## (-2.571)
## X_C2 -0.042*
## (-1.867)
## A3_K3 -0.024
## (-0.802)
## C3_A3 -0.009**
## (-2.564)
## X_C3 -0.037
## (-1.688)
## ---------------------------------------------------
## adj. R-squared 0.972 0.979 0.968
## N 27 27 27
## ===================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(akiba1_X_ols)##
## Durbin-Watson test
##
## data: akiba1_X_ols
## DW = 1.7755, p-value = 0.02002
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba2_X_ols)##
## Durbin-Watson test
##
## data: akiba2_X_ols
## DW = 1.687, p-value = 0.01055
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba3_X_ols)##
## Durbin-Watson test
##
## data: akiba3_X_ols
## DW = 1.6088, p-value = 0.006017
## alternative hypothesis: true autocorrelation is greater than 0akiba1_XG2_ols <- lm(K1_N~A1_K1+C1_A1+X_C1+U_L+W+G2_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba2_XG2_ols <- lm(K2_N~A2_K2+C2_A2+X_C2+U_L+W+G2_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba3_XG2_ols <- lm(K3_N~A3_K3+C3_A3+X_C3+U_L+W+G2_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba_XG2_model <- mtable("1"=akiba1_XG2_ols,"2"=akiba2_XG2_ols,"3"=akiba3_XG2_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_XG2_model)##
## Calls:
## 1: lm(formula = K1_N ~ A1_K1 + C1_A1 + X_C1 + U_L + W + G2_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + U_L + W + G2_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 3: lm(formula = K3_N ~ A3_K3 + C3_A3 + X_C3 + U_L + W + G2_N + F1_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
##
## ===================================================
## 1 2 3
## ---------- ---------- ----------
## K1_N K2_N K3_N
## ---------------------------------------------------
## (Intercept) 3.769 3.079 4.936
## (1.145) (0.852) (1.246)
## A1_K1 -0.032
## (-1.122)
## C1_A1 -0.009**
## (-2.549)
## X_C1 -0.049*
## (-1.904)
## U_L 0.347** 0.346** 0.311**
## (2.872) (2.884) (2.778)
## W 0.006 0.005 0.005
## (0.609) (0.449) (0.573)
## G2_N -0.002 -0.003 -0.004
## (-0.358) (-0.483) (-0.676)
## F1_N 0.097 0.132 0.085
## (0.970) (1.318) (0.931)
## M1524_N 0.037 -0.002 -0.002
## (0.282) (-0.011) (-0.017)
## E2_N 0.046 0.126 -0.106
## (0.128) (0.345) (-0.290)
## Y1972 0.100 0.110 0.086
## (0.823) (0.877) (0.734)
## TREND -0.021 -0.012 0.002
## (-0.216) (-0.118) (0.024)
## A2_K2 -0.024
## (-0.779)
## C2_A2 -0.010**
## (-2.642)
## X_C2 -0.047*
## (-2.007)
## A3_K3 -0.037
## (-1.160)
## C3_A3 -0.009**
## (-2.696)
## X_C3 -0.045*
## (-1.887)
## ---------------------------------------------------
## adj. R-squared 0.972 0.979 0.969
## N 27 27 27
## ===================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(akiba1_XG2_ols)##
## Durbin-Watson test
##
## data: akiba1_XG2_ols
## DW = 1.749, p-value = 0.01754
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba2_XG2_ols)##
## Durbin-Watson test
##
## data: akiba2_XG2_ols
## DW = 1.5814, p-value = 0.003929
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba3_XG2_ols)##
## Durbin-Watson test
##
## data: akiba3_XG2_ols
## DW = 1.5922, p-value = 0.004682
## alternative hypothesis: true autocorrelation is greater than 0akiba_XF2_ols <- lm(K2_N~A2_K2+C2_A2+X_C2+U_L+W+G1_N+F2_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba_XF3_ols <- lm(K2_N~A2_K2+C2_A2+X_C2+U_L+W+G1_N+F3_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
# 変数F3は1955年からデータを得ることができる
akiba5586 <- subset(capital, 1955<=YEAR & YEAR<=1986)
akiba5586$TREND=1:nrow(akiba5586) # トレンド変数
akiba_XF3_ols5586 <- lm(K2_N~A2_K2+C2_A2+X_C2+U_L+W+G1_N+F3_N+M1524_N+E2_N+Y1972+TREND, data=akiba5586)
akiba_XF23_model <- mtable("1"=akiba_XF2_ols,"2"=akiba_XF3_ols,"3"=akiba_XF3_ols5586,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_XF23_model)##
## Calls:
## 1: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + U_L + W + G1_N + F2_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + U_L + W + G1_N + F3_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba)
## 3: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + U_L + W + G1_N + F3_N +
## M1524_N + E2_N + Y1972 + TREND, data = akiba5586)
##
## ====================================================
## 1 2 3
## ----------------------------------------------------
## (Intercept) -1.296 0.389 5.423
## (-0.308) (0.090) (0.756)
## A2_K2 -0.012 -0.007 -0.001
## (-0.413) (-0.231) (-0.013)
## C2_A2 -0.015*** -0.010** -0.009
## (-3.482) (-2.562) (-1.350)
## X_C2 -0.038* -0.035 -0.003
## (-1.801) (-1.567) (-0.109)
## U_L 0.390*** 0.355** 0.265**
## (3.840) (2.836) (2.233)
## W 0.015 0.023 -0.016
## (1.042) (1.273) (-0.566)
## G1_N -0.000 -0.003 0.004
## (-0.102) (-0.444) (0.450)
## F2_N 0.309*
## (1.923)
## M1524_N 0.018 0.127* -0.023
## (0.284) (1.771) (-0.206)
## E2_N 0.806 0.056 -0.412
## (1.466) (0.222) (-1.101)
## Y1972 0.125 0.064 0.114
## (1.001) (0.512) (0.596)
## TREND -0.026 -0.048 -0.093
## (-0.565) (-0.978) (-1.392)
## F3_N 0.116 -0.098
## (1.309) (-0.821)
## ----------------------------------------------------
## adj. R-squared 0.981 0.979 0.964
## N 27 27 32
## ====================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(akiba_XF2_ols)##
## Durbin-Watson test
##
## data: akiba_XF2_ols
## DW = 2.0271, p-value = 0.05079
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba_XF3_ols)##
## Durbin-Watson test
##
## data: akiba_XF3_ols
## DW = 1.7799, p-value = 0.009729
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba_XF3_ols5586)##
## Durbin-Watson test
##
## data: akiba_XF3_ols5586
## DW = 1.4307, p-value = 0.0007794
## alternative hypothesis: true autocorrelation is greater than 0変数の変換
# ラグ
akiba$K2LAG <- c(NA,akiba$K2[1:length(akiba$K2)-1])
akiba$A2_K2LAG <- (akiba$A2)/(akiba$K2LAG)
akiba$LA2_K2LAG <- log(akiba$A2_K2LAG)
# 対数
akiba$LK2_N <- log(akiba$K2_N)
akiba$LA2_K2 <- log(akiba$A2_K2)
akiba$LC2_A2 <- log(akiba$C2_A2)
akiba$LU_L <- log(akiba$U_L)
akiba$LW <- log(akiba$W)
akiba$LG1_N <- log(akiba$G1_N)
akiba$LF1_N <- log(akiba$F1_N)
akiba$LM1524_N <- log(akiba$M1524_N)
akiba$LE2_N <- log(akiba$E2_N)
#対数変換用の0の処理1(X=0をX=1に、Ehrlich 1975の方法)
akiba$X_C2_1 <- ifelse(akiba$X==0,1,akiba$X)/akiba$C2*100
akiba$LX_C2_1 <- log(akiba$X_C2_1) #対数変換での0の処理2(Xすべてに1を足す)
akiba$X_C2_2 <- (akiba$X+1)/akiba$C2*100
akiba$LX_C2_2 <- log(akiba$X_C2_2) 表2-8のうち「(2)線形ラグ」
akiba_lag_ols <- lm(K2_N~A2_K2LAG+C2_A2+X_C2+U_L+W+G1_N+F1_N+M1524_N+E2_N+Y1972+TREND, data=akiba)
akiba_lag_model <- mtable("1"=akiba_lag_ols,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_lag_model)##
## Calls:
## 1: lm(formula = K2_N ~ A2_K2LAG + C2_A2 + X_C2 + U_L + W + G1_N +
## F1_N + M1524_N + E2_N + Y1972 + TREND, data = akiba)
##
## =============================
## (Intercept) -1.403
## (-0.856)
## A2_K2LAG 0.311
## (0.844)
## C2_A2 -0.008
## (-1.549)
## X_C2 -0.031
## (-1.445)
## U_L 0.310**
## (2.246)
## W 0.016
## (1.335)
## G1_N -0.004
## (-0.906)
## F1_N 0.190*
## (1.948)
## M1524_N 0.064
## (1.164)
## E2_N 0.512
## (1.297)
## Y1972 0.088
## (0.704)
## TREND -0.015
## (-0.321)
## -----------------------------
## adj. R-squared 0.977
## N 26
## =============================
## Significance:
## *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1表2-8のうち「(1)対数」と「(3)対数ラグ」
akiba_log_ols1 <- lm(LK2_N~LA2_K2+LC2_A2+LX_C2_1+LU_L+LW+LG1_N+LF1_N+LM1524_N+LE2_N+Y1972+TREND, data=akiba)
akiba_laglog_ols1 <- lm(LK2_N~LA2_K2LAG+LC2_A2+LX_C2_1+LU_L+LW+LG1_N+LF1_N+LM1524_N+LE2_N+Y1972+TREND, data=akiba)
akiba_log_model1 <- mtable("1"=akiba_log_ols1,"2"=akiba_laglog_ols1,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_log_model1)##
## Calls:
## 1: lm(formula = LK2_N ~ LA2_K2 + LC2_A2 + LX_C2_1 + LU_L + LW +
## LG1_N + LF1_N + LM1524_N + LE2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = LK2_N ~ LA2_K2LAG + LC2_A2 + LX_C2_1 + LU_L + LW +
## LG1_N + LF1_N + LM1524_N + LE2_N + Y1972 + TREND, data = akiba)
##
## ======================================
## 1 2
## --------------------------------------
## (Intercept) -2.506 1.729
## (-0.285) (0.959)
## LA2_K2 0.924
## (0.454)
## LC2_A2 -0.214 -0.130
## (-1.649) (-0.798)
## LX_C2_1 -0.004 -0.002
## (-0.401) (-0.165)
## LU_L 0.224 0.190
## (1.579) (1.231)
## LW 0.462 0.536
## (1.039) (1.491)
## LG1_N -0.617 -0.773
## (-0.826) (-1.357)
## LF1_N -0.115 0.130
## (-0.198) (0.230)
## LM1524_N 0.704 0.474
## (1.733) (1.310)
## LE2_N 0.078 0.165
## (0.230) (0.563)
## Y1972 0.020 0.002
## (0.264) (0.030)
## TREND -0.021 -0.010
## (-0.717) (-0.402)
## LA2_K2LAG 0.281
## (1.269)
## --------------------------------------
## adj. R-squared 0.969 0.970
## N 27 26
## ======================================
## Significance: *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(akiba_lag_ols)##
## Durbin-Watson test
##
## data: akiba_lag_ols
## DW = 2.068, p-value = 0.09291
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba_log_ols1)##
## Durbin-Watson test
##
## data: akiba_log_ols1
## DW = 1.8087, p-value = 0.03568
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba_laglog_ols1)##
## Durbin-Watson test
##
## data: akiba_laglog_ols1
## DW = 1.7299, p-value = 0.01866
## alternative hypothesis: true autocorrelation is greater than 0akiba_log_ols2 <- lm(LK2_N~LA2_K2+LC2_A2+LX_C2_2+LU_L+LW+LG1_N+LF1_N+LM1524_N+LE2_N+Y1972+TREND, data=akiba)
akiba_laglog_ols2 <- lm(LK2_N~LA2_K2LAG+LC2_A2+LX_C2_2+LU_L+LW+LG1_N+LF1_N+LM1524_N+LE2_N+Y1972+TREND, data=akiba)
akiba_log_model2 <- mtable("1"=akiba_log_ols2,"2"=akiba_laglog_ols2,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_log_model2)##
## Calls:
## 1: lm(formula = LK2_N ~ LA2_K2 + LC2_A2 + LX_C2_2 + LU_L + LW +
## LG1_N + LF1_N + LM1524_N + LE2_N + Y1972 + TREND, data = akiba)
## 2: lm(formula = LK2_N ~ LA2_K2LAG + LC2_A2 + LX_C2_2 + LU_L + LW +
## LG1_N + LF1_N + LM1524_N + LE2_N + Y1972 + TREND, data = akiba)
##
## ======================================
## 1 2
## --------------------------------------
## (Intercept) -2.014 1.526
## (-0.233) (0.852)
## LA2_K2 0.763
## (0.381)
## LC2_A2 -0.212 -0.142
## (-1.661) (-0.874)
## LX_C2_2 -0.008 -0.005
## (-0.804) (-0.497)
## LU_L 0.239 0.211
## (1.716) (1.372)
## LW 0.475 0.568
## (1.091) (1.632)
## LG1_N -0.593 -0.774
## (-0.805) (-1.376)
## LF1_N -0.060 0.171
## (-0.106) (0.309)
## LM1524_N 0.675 0.486
## (1.681) (1.362)
## LE2_N 0.083 0.172
## (0.251) (0.599)
## Y1972 0.016 -0.000
## (0.219) (-0.003)
## TREND -0.021 -0.011
## (-0.743) (-0.415)
## LA2_K2LAG 0.251
## (1.125)
## --------------------------------------
## adj. R-squared 0.970 0.970
## N 27 26
## ======================================
## Significance: *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(akiba_log_ols2)##
## Durbin-Watson test
##
## data: akiba_log_ols2
## DW = 1.7491, p-value = 0.02725
## alternative hypothesis: true autocorrelation is greater than 0dwtest(akiba_laglog_ols2)##
## Durbin-Watson test
##
## data: akiba_laglog_ols2
## DW = 1.7406, p-value = 0.02075
## alternative hypothesis: true autocorrelation is greater than 0注65：分散不均一性の検定 (p.348)
bptest(akiba1_X_ols)##
## studentized Breusch-Pagan test
##
## data: akiba1_X_ols
## BP = 12.898, df = 11, p-value = 0.3001bptest(akiba2_X_ols)##
## studentized Breusch-Pagan test
##
## data: akiba2_X_ols
## BP = 9.4742, df = 11, p-value = 0.5782bptest(akiba3_X_ols)##
## studentized Breusch-Pagan test
##
## data: akiba3_X_ols
## BP = 11.682, df = 11, p-value = 0.388akiba_cor <- cor(akiba[,c("A2_K2","C2_A2","X_C2","U_L","W","G1_N","F1_N","M1524_N","E2_N","Y1972","TREND")])
print(round(akiba_cor,3))## A2_K2 C2_A2 X_C2 U_L W G1_N F1_N M1524_N E2_N Y1972
## A2_K2 1.000 0.054 -0.323 0.190 -0.058 -0.057 0.088 -0.162 -0.207 0.071
## C2_A2 0.054 1.000 -0.178 0.107 -0.148 -0.117 0.155 0.089 -0.266 -0.088
## X_C2 -0.323 -0.178 1.000 -0.414 -0.431 -0.445 0.466 0.426 -0.348 -0.492
## U_L 0.190 0.107 -0.414 1.000 0.783 0.801 -0.750 -0.879 0.600 0.756
## W -0.058 -0.148 -0.431 0.783 1.000 0.984 -0.983 -0.913 0.943 0.930
## G1_N -0.057 -0.117 -0.445 0.801 0.984 1.000 -0.983 -0.871 0.930 0.870
## F1_N 0.088 0.155 0.466 -0.750 -0.983 -0.983 1.000 0.868 -0.968 -0.888
## M1524_N -0.162 0.089 0.426 -0.879 -0.913 -0.871 0.868 1.000 -0.780 -0.917
## E2_N -0.207 -0.266 -0.348 0.600 0.943 0.930 -0.968 -0.780 1.000 0.829
## Y1972 0.071 -0.088 -0.492 0.756 0.930 0.870 -0.888 -0.917 0.829 1.000
## TREND -0.014 -0.063 -0.482 0.861 0.967 0.991 -0.971 -0.887 0.892 0.861
## TREND
## A2_K2 -0.014
## C2_A2 -0.063
## X_C2 -0.482
## U_L 0.861
## W 0.967
## G1_N 0.991
## F1_N -0.971
## M1524_N -0.887
## E2_N 0.892
## Y1972 0.861
## TREND 1.000library(car)
akiba2_X_vif <- vif(akiba2_X_ols)
print(round(akiba2_X_vif,3))## A2_K2 C2_A2 X_C2 U_L W G1_N F1_N M1524_N
## 2.614 1.474 2.622 38.200 540.472 1045.290 591.224 40.888
## E2_N Y1972 TREND
## 156.874 30.646 1152.965akiba_step <- lm(K2_N~A2_K2+C2_A2+X_C2+F1_N+M1524_N, data=akiba)
akiba_step_model <- mtable("1"=akiba_step,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(akiba_step_model)##
## Calls:
## 1: lm(formula = K2_N ~ A2_K2 + C2_A2 + X_C2 + F1_N + M1524_N, data = akiba)
##
## ==============================
## (Intercept) 6.197*
## (1.743)
## A2_K2 -0.045
## (-1.261)
## C2_A2 -0.012**
## (-2.444)
## X_C2 -0.038
## (-1.435)
## F1_N 0.156***
## (11.825)
## M1524_N -0.076**
## (-2.450)
## ------------------------------
## adj. R-squared 0.958
## N 27
## ==============================
## Significance:
## *** = p < 0.01;
## ** = p < 0.05;
## * = p < 0.1ダービンワトソン比
dwtest(akiba_step)##
## Durbin-Watson test
##
## data: akiba_step
## DW = 1.0566, p-value = 0.0004979
## alternative hypothesis: true autocorrelation is greater than 0akiba_step_vif <- vif(akiba_step)
print(akiba_step_vif)## A2_K2 C2_A2 X_C2 F1_N M1524_N
## 1.678158 1.151086 1.779411 6.695957 5.592502トレンドと定数項
library(CADFtest)
CADFtest(akiba$K2_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$K2_N
## ADF(0) = -3.0274, p-value = 0.1474
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.563594CADFtest(akiba$A2_K2, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$A2_K2
## ADF(3) = -2.1759, p-value = 0.4787
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.8246995CADFtest(akiba$C2_A2, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$C2_A2
## ADF(3) = -1.7019, p-value = 0.7159
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.8069091CADFtest(akiba$X_C2, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$X_C2
## ADF(0) = -3.7814, p-value = 0.0376
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.7866579CADFtest(akiba$U_L, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$U_L
## ADF(4) = -1.9424, p-value = 0.5989
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.3045104CADFtest(akiba$W, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$W
## ADF(1) = -1.1674, p-value = 0.8925
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.08290479CADFtest(akiba$G1_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$G1_N
## ADF(0) = -1.5242, p-value = 0.7894
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1668875CADFtest(akiba$F1_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$F1_N
## ADF(0) = -0.46608, p-value = 0.9773
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.03379476CADFtest(akiba$M1524_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$M1524_N
## ADF(2) = -0.73489, p-value = 0.957
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.05088188CADFtest(akiba$E2_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$E2_N
## ADF(4) = -1.8896, p-value = 0.6257
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.08050493定数項
CADFtest(akiba$K2_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$K2_N
## ADF(0) = -1.6687, p-value = 0.4324
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.09503703CADFtest(akiba$A2_K2, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$A2_K2
## ADF(3) = -1.8449, p-value = 0.3504
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.6812186CADFtest(akiba$C2_A2, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$C2_A2
## ADF(1) = -2.9142, p-value = 0.05977
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.9629817CADFtest(akiba$X_C2, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$X_C2
## ADF(4) = -1.1026, p-value = 0.6958
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.4354885CADFtest(akiba$U_L, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$U_L
## ADF(0) = 0.071658, p-value = 0.9558
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.004465974CADFtest(akiba$W, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$W
## ADF(1) = -1.8405, p-value = 0.3523
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.03983275CADFtest(akiba$G1_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$G1_N
## ADF(1) = -1.372, p-value = 0.5769
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.02704914CADFtest(akiba$F1_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$F1_N
## ADF(0) = -4.3043, p-value = 0.003046
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.07686109CADFtest(akiba$M1524_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$M1524_N
## ADF(1) = -2.461, p-value = 0.138
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.05438632CADFtest(akiba$E2_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$E2_N
## ADF(1) = -4.2187, p-value = 0.00369
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.09532783どちらもなし
CADFtest(akiba$K2_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$K2_N
## ADF(0) = -3.0564, p-value = 0.003947
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.02808668CADFtest(akiba$A2_K2, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$A2_K2
## ADF(1) = -0.2956, p-value = 0.5676
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0004685171CADFtest(akiba$C2_A2, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$C2_A2
## ADF(1) = -0.034915, p-value = 0.6602
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0005616882CADFtest(akiba$X_C2, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$X_C2
## ADF(3) = -0.97943, p-value = 0.2828
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1951497CADFtest(akiba$U_L, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$U_L
## ADF(0) = 2.1233, p-value = 0.9892
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.0389454CADFtest(akiba$W, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$W
## ADF(3) = 0.27325, p-value = 0.756
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.002344583CADFtest(akiba$G1_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$G1_N
## ADF(3) = 0.86376, p-value = 0.8893
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.009163559CADFtest(akiba$F1_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$F1_N
## ADF(3) = -2.0526, p-value = 0.04089
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.03969107CADFtest(akiba$M1524_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$M1524_N
## ADF(1) = -2.1553, p-value = 0.03278
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.008984078CADFtest(akiba$E2_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: akiba$E2_N
## ADF(1) = 0.24849, p-value = 0.749
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.001009297トレンドと定数項
library(urca)
summary(ur.kpss(akiba$K2_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1819
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$A2_K2,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.144
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$C2_A2,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1418
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$X_C2,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.0702
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$U_L,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2169
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$W,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1881
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$G1_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1965
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$F1_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2525
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$M1524_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1306
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(akiba$E2_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2497
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216定数項
summary(ur.kpss(akiba$K2_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9662
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$A2_K2,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.1426
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$C2_A2,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.1533
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$X_C2,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.4849
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$U_L,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.8236
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$W,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9501
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$G1_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9838
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$F1_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9656
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$M1524_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.8433
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(akiba$E2_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.8474
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739変数の作成
capital$Cm_K2 <- capital$Cm/capital$K2
capital$X_Cm <- capital$X/capital$Cm
capital$Gm_N <- capital$Gm/capital$N*100
capital$L_N <- capital$L/(capital$N*10)*100
capital$N2029_N <- capital$N2029/capital$N
capital$V_N <- capital$V/capital$N*1000
capital$N1524_N <- capital$N1524/capital$Nデータの年代をMerriman(1988)にそろえる
mer <- subset(capital, 1957<=YEAR & YEAR<=1982)変数の対数変換
mer$LK2_N <- log(mer$K2_N)
mer$LCm_K2 <- log(mer$Cm_K2)
mer$LU_L <- log(mer$U_L)
mer$LGm_N <- log(mer$Gm_N)
mer$LL_N <- log(mer$L_N)
mer$LN2029_N <- log(mer$N2029_N)
mer$LM2029_N <- log(mer$M2029_N)
#対数変換用の0の処理1(X=0をX=1に)
mer$X_Cm_1 <- ifelse(mer$X==0,1,mer$X)/mer$Cm
mer$LX_Cm_1 <- log(mer$X_Cm_1)その他の変数の作成
mer$TREND=1:nrow(mer)
mer$Gm_N_LAG <- c(NA,mer$Gm_N[1:length(mer$Gm_N)-1])
mer$DGm_N <- (mer$Gm_N-mer$Gm_N_LAG)/mer$Gm_N_LAG #変化率
mer$LDGm_N <- mer$LGm_N-log(mer$Gm_N_LAG) #対数変化率library(RcmdrMisc)
mer_sum <- numSummary(mer[,c("K2_N","Cm_K2","X_Cm","U_L","Gm_N","L_N","N2029_N"),drop=F],statistics=c("mean", "sd", "quantiles"), quantiles=c(0,.5,1))
mer_sum1 <- cbind(mer_sum$table,mer_sum$n)
print(round(mer_sum1,3))## mean sd 0% 50% 100%
## K2_N 2.107 0.461 1.439 1.982 2.916 26
## Cm_K2 0.565 0.043 0.478 0.566 0.675 26
## X_Cm 0.009 0.009 0.000 0.004 0.029 26
## U_L 1.626 0.428 1.118 1.425 2.355 26
## Gm_N 161.891 66.829 59.985 173.550 260.334 26
## L_N 0.485 0.007 0.476 0.484 0.499 26
## N2029_N 0.174 0.015 0.136 0.178 0.190 26OLS
mer_log1_ols1 <- lm(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)
mer_log1_ols2 <- lm(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+LN2029_N, data = mer)
mer_log1_ols3 <- lm(LK2_N~LCm_K2+LX_Cm_1+LU_L+LDGm_N+LL_N+TREND+LN2029_N, data = mer)
library(memisc)
mer_log1ols_model <- mtable("1"=mer_log1_ols1,"2"=mer_log1_ols2,"3"=mer_log1_ols3,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(mer_log1ols_model)##
## Calls:
## 1: lm(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N +
## TREND + LN2029_N, data = mer)
## 2: lm(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N +
## LN2029_N, data = mer)
## 3: lm(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LDGm_N + LL_N +
## TREND + LN2029_N, data = mer)
##
## ======================================================
## 1 2 3
## ------------------------------------------------------
## (Intercept) 0.184 3.822*** -0.395
## (0.051) (4.160) (-0.542)
## LCm_K2 -0.048 -0.032 -0.208*
## (-0.346) (-0.236) (-1.943)
## LX_Cm_1 -0.014 -0.015 -0.009
## (-1.486) (-1.560) (-1.403)
## LU_L 0.109 0.049 0.058
## (1.095) (0.594) (1.015)
## LGm_N 0.102 -0.397***
## (0.215) (-16.742)
## LL_N -0.974 0.010 -2.207**
## (-0.708) (0.010) (-2.722)
## TREND -0.035 -0.028***
## (-1.054) (-16.565)
## LN2029_N 0.199 0.700*** 0.196
## (0.386) (3.537) (1.348)
## LDGm_N 0.781***
## (2.941)
## ------------------------------------------------------
## adj. R-squared 0.966 0.966 0.983
## N 26 26 25
## ======================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(mer_log1_ols1)##
## Durbin-Watson test
##
## data: mer_log1_ols1
## DW = 1.2149, p-value = 0.0006123
## alternative hypothesis: true autocorrelation is greater than 0dwtest(mer_log1_ols2)##
## Durbin-Watson test
##
## data: mer_log1_ols2
## DW = 1.3519, p-value = 0.003367
## alternative hypothesis: true autocorrelation is greater than 0dwtest(mer_log1_ols3)##
## Durbin-Watson test
##
## data: mer_log1_ols3
## DW = 1.3629, p-value = 0.003011
## alternative hypothesis: true autocorrelation is greater than 0Prais-Winsten法
library(prais)
mer_log1_prais1 <- prais_winsten(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)
mer_log1_prais2 <- prais_winsten(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+LN2029_N, data = mer)
mer_log1_prais3 <- prais_winsten(LK2_N~LCm_K2+LX_Cm_1+LU_L+LDGm_N+LL_N+TREND+LN2029_N, data = mer)Prais-Winsten法での(1)のモデル
summary(mer_log1_prais1)##
## Call:
## prais_winsten(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.071401 -0.023727 0.009096 0.024390 0.066940
##
## AR(1) coefficient rho after 19 Iterations: 0.461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.445409 3.238699 -0.138 0.8921
## LCm_K2 -0.147699 0.112207 -1.316 0.2046
## LX_Cm_1 -0.016904 0.007193 -2.350 0.0304 *
## LU_L 0.161222 0.098492 1.637 0.1190
## LGm_N 0.238650 0.422558 0.565 0.5792
## LL_N -0.816586 1.425087 -0.573 0.5737
## TREND -0.044445 0.029715 -1.496 0.1521
## LN2029_N 0.148879 0.487281 0.306 0.7635
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03655 on 18 degrees of freedom
## Multiple R-squared: 0.9613, Adjusted R-squared: 0.9463
## F-statistic: 63.9 on 7 and 18 DF, p-value: 2.014e-11
##
## Durbin-Watson statistic (original): 1.215
## Durbin-Watson statistic (transformed): 1.463Prais-Winsten法での(2)のモデル
summary(mer_log1_prais2)##
## Call:
## prais_winsten(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08077 -0.01980 0.00768 0.01941 0.07410
##
## AR(1) coefficient rho after 20 Iterations: 0.3597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.038922 1.022543 3.950 0.000859 ***
## LCm_K2 -0.125769 0.120260 -1.046 0.308772
## LX_Cm_1 -0.017393 0.007828 -2.222 0.038633 *
## LU_L 0.075663 0.087912 0.861 0.400152
## LGm_N -0.395499 0.029762 -13.289 4.54e-11 ***
## LL_N 0.279917 1.158227 0.242 0.811619
## LN2029_N 0.763891 0.216576 3.527 0.002252 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.038 on 19 degrees of freedom
## Multiple R-squared: 0.9618, Adjusted R-squared: 0.9497
## F-statistic: 79.66 on 6 and 19 DF, p-value: 1.921e-12
##
## Durbin-Watson statistic (original): 1.352
## Durbin-Watson statistic (transformed): 1.549Prais-Winsten法での(3)のモデル
summary(mer_log1_prais3)##
## Call:
## prais_winsten(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LDGm_N +
## LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047997 -0.014791 -0.001091 0.014248 0.046265
##
## AR(1) coefficient rho after 10 Iterations: 0.3231
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.101579 0.777864 -0.131 0.8976
## LCm_K2 -0.196721 0.093693 -2.100 0.0510 .
## LX_Cm_1 -0.010774 0.005757 -1.871 0.0786 .
## LU_L 0.076792 0.061695 1.245 0.2301
## LDGm_N 0.736844 0.244411 3.015 0.0078 **
## LL_N -1.847362 0.883264 -2.092 0.0518 .
## TREND -0.027444 0.001812 -15.146 2.65e-11 ***
## LN2029_N 0.219515 0.159862 1.373 0.1875
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02651 on 17 degrees of freedom
## Multiple R-squared: 0.9838, Adjusted R-squared: 0.9772
## F-statistic: 147.8 on 7 and 17 DF, p-value: 5.745e-14
##
## Durbin-Watson statistic (original): 1.363
## Durbin-Watson statistic (transformed): 1.62OLS
mer_ols1 <- lm(K2_N~Cm_K2+X_Cm+U_L+Gm_N+L_N+TREND+N2029_N, data = mer)
mer_ols2 <- lm(K2_N~Cm_K2+X_Cm+U_L+Gm_N+L_N+N2029_N, data = mer)
mer_ols3 <- lm(K2_N~Cm_K2+X_Cm+U_L+DGm_N+L_N+TREND+N2029_N, data = mer)
library(memisc)
mer_ols_model <- mtable("1"=mer_ols1,"2"=mer_ols2,"3"=mer_ols3,coef.style="stat",signif.symbols=c("***"=.01,"**"=.05, "*"=.1),summary.stats=c("adj. R-squared","N"))
print(mer_ols_model)##
## Calls:
## 1: lm(formula = K2_N ~ Cm_K2 + X_Cm + U_L + Gm_N + L_N + TREND +
## N2029_N, data = mer)
## 2: lm(formula = K2_N ~ Cm_K2 + X_Cm + U_L + Gm_N + L_N + N2029_N,
## data = mer)
## 3: lm(formula = K2_N ~ Cm_K2 + X_Cm + U_L + DGm_N + L_N + TREND +
## N2029_N, data = mer)
##
## =======================================================
## 1 2 3
## -------------------------------------------------------
## (Intercept) 3.902 3.235 7.951***
## (1.723) (1.474) (4.896)
## Cm_K2 -0.744 -0.887* -1.089***
## (-1.467) (-1.800) (-3.087)
## X_Cm -4.286* -4.143* -2.865*
## (-1.892) (-1.821) (-1.793)
## U_L 0.192** 0.214** 0.125*
## (2.128) (2.424) (2.093)
## Gm_N -0.004 -0.007***
## (-1.224) (-19.104)
## L_N -2.418 -1.421 -9.518***
## (-0.635) (-0.382) (-3.330)
## TREND -0.031 -0.064***
## (-1.101) (-20.410)
## N2029_N 3.166 5.184** -0.035
## (1.139) (2.470) (-0.023)
## DGm_N 1.471***
## (3.123)
## -------------------------------------------------------
## adj. R-squared 0.973 0.972 0.987
## N 26 26 25
## =======================================================
## Significance: *** = p < 0.01; ** = p < 0.05;
## * = p < 0.1ダービン・ワトソン比
dwtest(mer_ols1)##
## Durbin-Watson test
##
## data: mer_ols1
## DW = 1.2052, p-value = 0.0002535
## alternative hypothesis: true autocorrelation is greater than 0dwtest(mer_ols2)##
## Durbin-Watson test
##
## data: mer_ols2
## DW = 1.2254, p-value = 0.001074
## alternative hypothesis: true autocorrelation is greater than 0dwtest(mer_ols3)##
## Durbin-Watson test
##
## data: mer_ols3
## DW = 1.609, p-value = 0.02135
## alternative hypothesis: true autocorrelation is greater than 0Prais-Winsten法
mer_prais1 <- prais_winsten(K2_N~Cm_K2+X_Cm+U_L+Gm_N+L_N+TREND+N2029_N, data = mer)
mer_prais2 <- prais_winsten(K2_N~Cm_K2+X_Cm+U_L+Gm_N+L_N+N2029_N, data = mer)
mer_prais3 <- prais_winsten(K2_N~Cm_K2+X_Cm+U_L+DGm_N+L_N+TREND+N2029_N, data = mer)Prais-Winsten法での(1)のモデル
summary(mer_prais1)##
## Call:
## prais_winsten(formula = K2_N ~ Cm_K2 + X_Cm + U_L + Gm_N + L_N +
## TREND + N2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.123498 -0.035520 -0.001859 0.024738 0.188755
##
## AR(1) coefficient rho after 13 Iterations: 0.3967
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.274780 2.479454 1.321 0.2031
## Cm_K2 -0.931165 0.427209 -2.180 0.0428 *
## X_Cm -5.027813 1.785530 -2.816 0.0114 *
## U_L 0.205798 0.094865 2.169 0.0437 *
## Gm_N -0.003405 0.003520 -0.967 0.3461
## L_N -1.064368 4.379921 -0.243 0.8107
## TREND -0.033732 0.032272 -1.045 0.3097
## N2029_N 3.404385 3.134480 1.086 0.2918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0709 on 18 degrees of freedom
## Multiple R-squared: 0.9742, Adjusted R-squared: 0.9641
## F-statistic: 97.01 on 7 and 18 DF, p-value: 5.463e-13
##
## Durbin-Watson statistic (original): 1.205
## Durbin-Watson statistic (transformed): 1.508Prais-Winsten法での(2)のモデル
summary(mer_prais2)##
## Call:
## prais_winsten(formula = K2_N ~ Cm_K2 + X_Cm + U_L + Gm_N + L_N +
## N2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09398 -0.03153 -0.01453 0.03412 0.19364
##
## AR(1) coefficient rho after 13 Iterations: 0.4199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4230734 2.3880891 1.015 0.3230
## Cm_K2 -1.0660727 0.4074404 -2.617 0.0170 *
## X_Cm -4.9411178 1.7650477 -2.799 0.0114 *
## U_L 0.2264578 0.0936766 2.417 0.0258 *
## Gm_N -0.0070432 0.0004671 -15.078 5.03e-12 ***
## L_N 0.2852685 4.2763262 0.067 0.9475
## N2029_N 5.5520530 2.4397945 2.276 0.0346 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07099 on 19 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9632
## F-statistic: 110.2 on 6 and 19 DF, p-value: 9.918e-14
##
## Durbin-Watson statistic (original): 1.225
## Durbin-Watson statistic (transformed): 1.61Prais-Winsten法での(3)のモデル
summary(mer_prais3)##
## Call:
## prais_winsten(formula = K2_N ~ Cm_K2 + X_Cm + U_L + DGm_N + L_N +
## TREND + N2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.074510 -0.031614 -0.003724 0.026923 0.079876
##
## AR(1) coefficient rho after 16 Iterations: 0.2165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.339156 1.705540 4.303 0.000482 ***
## Cm_K2 -0.984568 0.328437 -2.998 0.008094 **
## X_Cm -3.252575 1.490678 -2.182 0.043439 *
## U_L 0.136458 0.063259 2.157 0.045601 *
## DGm_N 1.417693 0.453660 3.125 0.006165 **
## L_N -8.491412 3.066510 -2.769 0.013132 *
## TREND -0.063415 0.003315 -19.127 6.2e-13 ***
## N2029_N 0.173593 1.634103 0.106 0.916642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05055 on 17 degrees of freedom
## Multiple R-squared: 0.9892, Adjusted R-squared: 0.9848
## F-statistic: 223.5 on 7 and 17 DF, p-value: 1.818e-15
##
## Durbin-Watson statistic (original): 1.609
## Durbin-Watson statistic (transformed): 1.687Breusch-Godfrey検定 (p.368)
bgtest(mer_log1_ols1,order=2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: mer_log1_ols1
## LM test = 5.573, df = 2, p-value = 0.06164bgtest(mer_log1_ols2,order=2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: mer_log1_ols2
## LM test = 4.058, df = 2, p-value = 0.1315bgtest(mer_log1_ols3,order=2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: mer_log1_ols3
## LM test = 5.3587, df = 2, p-value = 0.06861mer_log1_cor <- cor(mer[,c("LCm_K2","LX_Cm_1","LU_L","LGm_N","LDGm_N","LL_N","TREND","LN2029_N")],use="pairwise.complete.obs")
print(round(mer_log1_cor,3))## LCm_K2 LX_Cm_1 LU_L LGm_N LDGm_N LL_N TREND LN2029_N
## LCm_K2 1.000 0.083 -0.288 -0.526 0.374 0.016 -0.531 0.259
## LX_Cm_1 0.083 1.000 -0.168 -0.420 0.244 -0.021 -0.481 0.601
## LU_L -0.288 -0.168 1.000 0.177 -0.479 -0.727 0.372 -0.677
## LGm_N -0.526 -0.420 0.177 1.000 -0.597 0.127 0.974 -0.491
## LDGm_N 0.374 0.244 -0.479 -0.597 1.000 0.472 -0.663 0.434
## LL_N 0.016 -0.021 -0.727 0.127 0.472 1.000 -0.023 0.225
## TREND -0.531 -0.481 0.372 0.974 -0.663 -0.023 1.000 -0.659
## LN2029_N 0.259 0.601 -0.677 -0.491 0.434 0.225 -0.659 1.000(1)のモデルのVIF
mer_log1ols1_vif<-vif(mer_log1_ols1)
print(round(mer_log1ols1_vif,3))## LCm_K2 LX_Cm_1 LU_L LGm_N LL_N TREND LN2029_N
## 1.667 2.028 10.697 800.671 6.303 1021.123 35.762(2)のモデルのVIF
mer_log1ols2_vif<-vif(mer_log1_ols2)
print(round(mer_log1ols2_vif,3))## LCm_K2 LX_Cm_1 LU_L LGm_N LL_N LN2029_N
## 1.649 2.018 7.131 1.996 3.402 5.255(3)のモデルのVIF
mer_log1ols3_vif<-vif(mer_log1_ols3)
print(round(mer_log1ols3_vif,3))## LCm_K2 LX_Cm_1 LU_L LDGm_N LL_N TREND LN2029_N
## 2.024 1.899 7.205 2.797 4.418 4.699 5.878変数の準備
mer$LK2_NLAG <- c(NA,mer$LK2_N[1:length(mer$LK2_N)-1])
mer$LCm_K2LAG <- c(NA,mer$LCm_K2[1:length(mer$LCm_K2)-1])
mer$LX_Cm_1LAG <- c(NA,mer$LX_Cm_1[1:length(mer$LX_Cm_1)-1])
mer$LU_LLAG <- c(NA,mer$LU_L[1:length(mer$LU_L)-1])
mer$LGm_NLAG <- c(NA,mer$LGm_N[1:length(mer$LGm_N)-1])
mer$LL_NLAG <- c(NA,mer$LL_N[1:length(mer$LL_N)-1])
mer$LN2029_NLAG <- c(NA,mer$LN2029_N[1:length(mer$LN2029_N)-1])
mer$LV_N <- log(mer$V_N)
mer$LN <- log(mer$N)
mer$LN1524_N <- log(mer$N1524_N)library(estimatr)
iv_estm1_1 <- iv_robust(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N | LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N+LK2_NLAG+LCm_K2LAG+LX_Cm_1LAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG+LV_N+LN+LN1524_N,se_type = "classical",data=mer,diagnostics = TRUE)
summary(iv_estm1_1)##
## Call:
## iv_robust(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + TREND + LN2029_N | LX_Cm_1 + LU_L + LGm_N + LL_N +
## TREND + LN2029_N + LK2_NLAG + LCm_K2LAG + LX_Cm_1LAG + LU_LLAG +
## LGm_NLAG + LL_NLAG + LN2029_NLAG + LV_N + LN + LN1524_N,
## data = mer, se_type = "classical", diagnostics = TRUE)
##
## Standard error type: classical
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 1.717798 3.291872 0.52183 0.6085 -5.22744 8.663040 17
## LCm_K2 -0.006202 0.203851 -0.03043 0.9761 -0.43629 0.423886 17
## LX_Cm_1 -0.011027 0.008722 -1.26417 0.2232 -0.02943 0.007376 17
## LU_L 0.113545 0.094769 1.19813 0.2473 -0.08640 0.313489 17
## LGm_N -0.091430 0.437916 -0.20878 0.8371 -1.01535 0.832491 17
## LL_N -0.459567 1.280878 -0.35879 0.7242 -3.16198 2.242850 17
## TREND -0.022344 0.030711 -0.72754 0.4768 -0.08714 0.042451 17
## LN2029_N 0.387197 0.474387 0.81621 0.4257 -0.61367 1.388065 17
##
## Multiple R-squared: 0.9794 , Adjusted R-squared: 0.9709
## F-statistic: 115.4 on 7 and 17 DF, p-value: 4.499e-13
##
## Diagnostics:
## numdf dendf value p.value
## Weak instruments 10 8 0.730 0.686
## Wu-Hausman 1 16 2.551 0.130
## Overidentifying 9 NA 13.183 0.155操作変数からLCm_K2LAGを除いた場合
iv_estm1_2 <- iv_robust(LK2_N~LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N | LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N+LK2_NLAG+LX_Cm_1LAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG+LV_N+LN+LN1524_N,se_type = "classical",data=mer,diagnostics = TRUE)
summary(iv_estm1_2)##
## Call:
## iv_robust(formula = LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + TREND + LN2029_N | LX_Cm_1 + LU_L + LGm_N + LL_N +
## TREND + LN2029_N + LK2_NLAG + LX_Cm_1LAG + LU_LLAG + LGm_NLAG +
## LL_NLAG + LN2029_NLAG + LV_N + LN + LN1524_N, data = mer,
## se_type = "classical", diagnostics = TRUE)
##
## Standard error type: classical
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 1.750772 3.322579 0.52693 0.6050 -5.25926 8.760801 17
## LCm_K2 0.008187 0.206670 0.03961 0.9689 -0.42785 0.444223 17
## LX_Cm_1 -0.011069 0.008803 -1.25740 0.2256 -0.02964 0.007504 17
## LU_L 0.115905 0.095698 1.21115 0.2424 -0.08600 0.317810 17
## LGm_N -0.091031 0.441960 -0.20597 0.8393 -1.02349 0.841424 17
## LL_N -0.433978 1.293180 -0.33559 0.7413 -3.16235 2.294393 17
## TREND -0.022257 0.030995 -0.71809 0.4825 -0.08765 0.043136 17
## LN2029_N 0.393327 0.478841 0.82141 0.4228 -0.61694 1.403593 17
##
## Multiple R-squared: 0.979 , Adjusted R-squared: 0.9704
## F-statistic: 113.3 on 7 and 17 DF, p-value: 5.238e-13
##
## Diagnostics:
## numdf dendf value p.value
## Weak instruments 9 9 0.897 0.563
## Wu-Hausman 1 16 2.920 0.107
## Overidentifying 8 NA 12.154 0.144操作変数からLCm_K2LAGを除いた場合
iv_estm2_2 <- iv_robust(LK2_N~LCm_K2+LU_L+LGm_N+LL_N+TREND+LN2029_N | LU_L+LGm_N+LL_N+TREND+LN2029_N+LK2_NLAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG+LV_N+LN+LN1524_N,se_type = "classical",data=mer,diagnostics = TRUE)
summary(iv_estm2_2)##
## Call:
## iv_robust(formula = LK2_N ~ LCm_K2 + LU_L + LGm_N + LL_N + TREND +
## LN2029_N | LU_L + LGm_N + LL_N + TREND + LN2029_N + LK2_NLAG +
## LU_LLAG + LGm_NLAG + LL_NLAG + LN2029_NLAG + LV_N + LN +
## LN1524_N, data = mer, se_type = "classical", diagnostics = TRUE)
##
## Standard error type: classical
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 1.08008 3.25161 0.3322 0.7436 -5.75131 7.91147 18
## LCm_K2 -0.04849 0.24590 -0.1972 0.8459 -0.56511 0.46814 18
## LU_L 0.07536 0.09278 0.8122 0.4273 -0.11957 0.27029 18
## LGm_N -0.06626 0.43458 -0.1525 0.8805 -0.97927 0.84675 18
## LL_N -0.84567 1.26576 -0.6681 0.5125 -3.50495 1.81360 18
## TREND -0.02433 0.03048 -0.7982 0.4352 -0.08838 0.03971 18
## LN2029_N 0.20966 0.45553 0.4603 0.6508 -0.74736 1.16669 18
##
## Multiple R-squared: 0.9785 , Adjusted R-squared: 0.9713
## F-statistic: 136.1 on 6 and 18 DF, p-value: 5.325e-14
##
## Diagnostics:
## numdf dendf value p.value
## Weak instruments 8 11 0.655 0.720
## Wu-Hausman 1 17 0.820 0.378
## Overidentifying 7 NA 11.957 0.102instr1_1 <- lm(LCm_K2~LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N +LK2_NLAG+LCm_K2LAG+LX_Cm_1LAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG +LV_N+LN+LN1524_N,data=mer)
resid1_1 <- c(NA,instr1_1$residuals)
endog1_1 <- lm(LK2_N~resid1_1+LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)
summary(endog1_1)##
## Call:
## lm(formula = LK2_N ~ resid1_1 + LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.051284 -0.011377 0.002231 0.020678 0.048874
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.717798 2.937572 0.585 0.567
## resid1_1 -0.401820 0.251559 -1.597 0.130
## LCm_K2 -0.006202 0.181911 -0.034 0.973
## LX_Cm_1 -0.011027 0.007784 -1.417 0.176
## LU_L 0.113545 0.084569 1.343 0.198
## LGm_N -0.091430 0.390783 -0.234 0.818
## LL_N -0.459567 1.143019 -0.402 0.693
## TREND -0.022344 0.027406 -0.815 0.427
## LN2029_N 0.387197 0.423329 0.915 0.374
##
## Residual standard error: 0.03245 on 16 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9846, Adjusted R-squared: 0.9768
## F-statistic: 127.5 on 8 and 16 DF, p-value: 5.132e-13library(lmtest)
WuHau1_1 <- waldtest(endog1_1, .~.-resid1_1)
print(WuHau1_1)## Wald test
##
## Model 1: LK2_N ~ resid1_1 + LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N + TREND +
## LN2029_N
## Model 2: LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N + TREND + LN2029_N
## Res.Df Df F Pr(>F)
## 1 16
## 2 17 -1 2.5514 0.1298Prais-Winsten法を使用
endog1_1_prais <- prais_winsten(LK2_N~resid1_1+LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)Prais-Winsten法の結果表示
summary(endog1_1_prais)##
## Call:
## prais_winsten(formula = LK2_N ~ resid1_1 + LCm_K2 + LX_Cm_1 +
## LU_L + LGm_N + LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.049880 -0.012569 0.001799 0.019935 0.051359
##
## AR(1) coefficient rho after 15 Iterations: 0.1242
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.181537 2.944018 0.401 0.693
## resid1_1 -0.356966 0.248229 -1.438 0.170
## LCm_K2 -0.034753 0.178656 -0.195 0.848
## LX_Cm_1 -0.011381 0.007498 -1.518 0.149
## LU_L 0.125585 0.086526 1.451 0.166
## LGm_N -0.013795 0.391506 -0.035 0.972
## LL_N -0.563051 1.172579 -0.480 0.638
## TREND -0.027894 0.027433 -1.017 0.324
## LN2029_N 0.315817 0.426646 0.740 0.470
##
## Residual standard error: 0.03227 on 16 degrees of freedom
## Multiple R-squared: 0.9822, Adjusted R-squared: 0.9733
## F-statistic: 110.4 on 8 and 16 DF, p-value: 1.586e-12
##
## Durbin-Watson statistic (original): 1.7
## Durbin-Watson statistic (transformed): 1.741操作変数LCm_K2LAGを除いた場合
instr1_2 <- lm(LCm_K2~LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N +LK2_NLAG+LX_Cm_1LAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG +LV_N+LN+LN1524_N,data=mer)
resid1_2 <- c(NA,instr1_2$residuals)
endog1_2 <- lm(LK2_N~resid1_2+LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)
summary(endog1_2)##
## Call:
## lm(formula = LK2_N ~ resid1_2 + LCm_K2 + LX_Cm_1 + LU_L + LGm_N +
## LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.053369 -0.010808 0.001773 0.021176 0.050688
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.750772 2.909107 0.602 0.556
## resid1_2 -0.425824 0.249207 -1.709 0.107
## LCm_K2 0.008187 0.180951 0.045 0.964
## LX_Cm_1 -0.011069 0.007708 -1.436 0.170
## LU_L 0.115905 0.083789 1.383 0.186
## LGm_N -0.091031 0.386961 -0.235 0.817
## LL_N -0.433978 1.132253 -0.383 0.707
## TREND -0.022257 0.027138 -0.820 0.424
## LN2029_N 0.393327 0.419253 0.938 0.362
##
## Residual standard error: 0.03213 on 16 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9849, Adjusted R-squared: 0.9773
## F-statistic: 130 on 8 and 16 DF, p-value: 4.389e-13WuHau1_2 <- waldtest(endog1_2, .~.-resid1_2)
print(WuHau1_2)## Wald test
##
## Model 1: LK2_N ~ resid1_2 + LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N + TREND +
## LN2029_N
## Model 2: LK2_N ~ LCm_K2 + LX_Cm_1 + LU_L + LGm_N + LL_N + TREND + LN2029_N
## Res.Df Df F Pr(>F)
## 1 16
## 2 17 -1 2.9197 0.1068Prais-Winsten法を使用
endog1_2_prais <- prais_winsten(LK2_N~resid1_2+LCm_K2+LX_Cm_1+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)Prais-Winsten法の結果表示
summary(endog1_2_prais)##
## Call:
## prais_winsten(formula = LK2_N ~ resid1_2 + LCm_K2 + LX_Cm_1 +
## LU_L + LGm_N + LL_N + TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.052100 -0.013242 0.001475 0.019038 0.052440
##
## AR(1) coefficient rho after 14 Iterations: 0.1011
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.307253 2.919086 0.448 0.660
## resid1_2 -0.388674 0.247539 -1.570 0.136
## LCm_K2 -0.014688 0.179190 -0.082 0.936
## LX_Cm_1 -0.011335 0.007489 -1.514 0.150
## LU_L 0.125319 0.085472 1.466 0.162
## LGm_N -0.027645 0.388233 -0.071 0.944
## LL_N -0.526413 1.156857 -0.455 0.655
## TREND -0.026789 0.027209 -0.985 0.339
## LN2029_N 0.334756 0.422564 0.792 0.440
##
## Residual standard error: 0.03201 on 16 degrees of freedom
## Multiple R-squared: 0.983, Adjusted R-squared: 0.9744
## F-statistic: 115.4 on 8 and 16 DF, p-value: 1.119e-12
##
## Durbin-Watson statistic (original): 1.741
## Durbin-Watson statistic (transformed): 1.772操作変数からLCm_K2LAGを除いた場合
instr2_2 <- lm(LCm_K2~LU_L+LGm_N+LL_N+TREND+LN2029_N +LK2_NLAG+LU_LLAG+LGm_NLAG+LL_NLAG+LN2029_NLAG +LV_N+LN+LN1524_N,data=mer)
resid2_2 <- c(NA,instr2_2$residuals)
endog2_2 <- lm(LK2_N~resid2_2+LCm_K2+LU_L+LGm_N+LL_N+TREND+LN2029_N, data = mer)
summary(endog2_2)##
## Call:
## lm(formula = LK2_N ~ resid2_2 + LCm_K2 + LU_L + LGm_N + LL_N +
## TREND + LN2029_N, data = mer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.051284 -0.013124 -0.001968 0.016589 0.054968
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.08008 3.12085 0.346 0.734
## resid2_2 -0.25970 0.28680 -0.906 0.378
## LCm_K2 -0.04849 0.23601 -0.205 0.840
## LU_L 0.07536 0.08905 0.846 0.409
## LGm_N -0.06626 0.41710 -0.159 0.876
## LL_N -0.84567 1.21486 -0.696 0.496
## TREND -0.02433 0.02926 -0.832 0.417
## LN2029_N 0.20966 0.43721 0.480 0.638
##
## Residual standard error: 0.03467 on 17 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9813, Adjusted R-squared: 0.9735
## F-statistic: 127.2 on 7 and 17 DF, p-value: 2.007e-13WuHau2_2 <- waldtest(endog2_2, .~.-resid2_2)
print(WuHau2_2)## Wald test
##
## Model 1: LK2_N ~ resid2_2 + LCm_K2 + LU_L + LGm_N + LL_N + TREND + LN2029_N
## Model 2: LK2_N ~ LCm_K2 + LU_L + LGm_N + LL_N + TREND + LN2029_N
## Res.Df Df F Pr(>F)
## 1 17
## 2 18 -1 0.8199 0.3779トレンドと定数項
CADFtest(mer$LK2_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LK2_N
## ADF(0) = -3.6819, p-value = 0.04666
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.7593317CADFtest(mer$LCm_K2, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LCm_K2
## ADF(1) = -2.801, p-value = 0.2119
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.9108214CADFtest(mer$LX_Cm_1, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LX_Cm_1
## ADF(1) = -2.4244, p-value = 0.3578
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.7592315CADFtest(mer$LU_L, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LU_L
## ADF(0) = -2.4265, p-value = 0.3568
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.3148128CADFtest(mer$LGm_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LGm_N
## ADF(0) = -0.55666, p-value = 0.9711
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0380436CADFtest(mer$LDGm_N, type="trend", max.lag.y=4, criterion="MAIC") ##
## ADF test
##
## data: mer$LDGm_N
## ADF(2) = -1.8533, p-value = 0.6404
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.704629CADFtest(mer$LL_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LL_N
## ADF(0) = -1.5455, p-value = 0.7797
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.2116962CADFtest(mer$LN2029_N, type="trend", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LN2029_N
## ADF(2) = -0.17394, p-value = 0.9891
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.01438287定数項
CADFtest(mer$LK2_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LK2_N
## ADF(2) = -1.2425, p-value = 0.6356
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.07782813CADFtest(mer$LCm_K2, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LCm_K2
## ADF(3) = -0.69682, p-value = 0.8268
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1970884CADFtest(mer$LX_Cm_1, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LX_Cm_1
## ADF(1) = -2.0727, p-value = 0.2565
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.6180107CADFtest(mer$LU_L, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LU_L
## ADF(0) = -0.020004, p-value = 0.9464
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.001866841CADFtest(mer$LGm_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LGm_N
## ADF(3) = -1.6775, p-value = 0.4275
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.04380674CADFtest(mer$LDGm_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LDGm_N
## ADF(2) = -1.0631, p-value = 0.709
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.2451331CADFtest(mer$LL_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LL_N
## ADF(0) = -1.3079, p-value = 0.606
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.161773CADFtest(mer$LN2029_N, type="drift", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LN2029_N
## ADF(2) = 0.31553, p-value = 0.9733
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.02679662どちらもなし
CADFtest(mer$LK2_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LK2_N
## ADF(0) = -3.0626, p-value = 0.003976
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.04693691CADFtest(mer$LCm_K2, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LCm_K2
## ADF(1) = 0.3169, p-value = 0.7678
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.008801401CADFtest(mer$LX_Cm_1, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LX_Cm_1
## ADF(1) = 0.1227, p-value = 0.7108
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.006446368CADFtest(mer$LU_L, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LU_L
## ADF(0) = 0.86978, p-value = 0.89
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.04205687CADFtest(mer$LGm_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LGm_N
## ADF(3) = 0.43477, p-value = 0.7989
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.001118678CADFtest(mer$LDGm_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LDGm_N
## ADF(2) = -1.3289, p-value = 0.1641
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.1341785CADFtest(mer$LL_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LL_N
## ADF(0) = -0.16694, p-value = 0.614
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## -0.0004205504CADFtest(mer$LN2029_N, type="none", max.lag.y=4, criterion="MAIC")##
## ADF test
##
## data: mer$LN2029_N
## ADF(1) = 1.0092, p-value = 0.9114
## alternative hypothesis: true delta is less than 0
## sample estimates:
## delta
## 0.002497968トレンドと定数項
library(urca)
summary(ur.kpss(mer$LK2_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.088
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LCm_K2,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.0884
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LX_Cm_1,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1017
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LU_L,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2218
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LGm_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2408
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LDGm_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.0933
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LL_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.1704
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216summary(ur.kpss(mer$LN2029_N,type="tau",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 0.2107
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216定数項
summary(ur.kpss(mer$LK2_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9291
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LCm_K2,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.5277
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LX_Cm_1,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.4088
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LU_L,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.3477
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LGm_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.9305
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LDGm_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.6243
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LL_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.1727
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739summary(ur.kpss(mer$LN2029_N,type="mu",lags="short"))##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: mu with 2 lags.
##
## Value of test-statistic is: 0.5232
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.347 0.463 0.574 0.739