パネルデータ分析のメモ。その２
【目次】
以下で基本部分を取りまとめた、今回は実データで推定してみる。
【統計】時系列分析（パネルデータ分析） - こちにぃるの日記
使用するデータ
Rのplmパッケージの Grunfeld を用いる。
目的変数はinv、説明はvalue captal、属性情報はfirm、経時情報はyearに含まれる。
・firm ：企業（ファーム）を表す変数。
・year ：年度。
・inv ：企業ごとの投資額。
・value ：企業の市場価値を表す変数。
・capital：企業の設備や設備の在庫を表す変数。
10社の1935年～1954年のデータが含まれる。
1950年までのデータを使ってパネルデータ分析をしてみる。
library(plm) data("Grunfeld", package="plm") Grunfeld$FLG <- ifelse(Grunfeld$year<=1950,1,0) ADTS_T <- pdata.frame(Grunfeld[Grunfeld$FLG==1,], index=c("firm","year"))
パネルデータ分析の実装
固定効果モデルの実行。
valueとcapitalの両方が統計的に有意である。
価値valueが1単位増加すると、投資 inv は約0.0752単位増加する
資本capitalが1単位増加すると、投資 inv は約0.2063単位増加する
> #-- 固定効果モデル > fit1 <- plm(inv ~ value + capital, data=ADTS_T, model="within") > summary(fit1) Oneway (individual) effect Within Model Call: plm(formula = inv ~ value + capital, data = ADTS_T, model = "within") Balanced Panel: n = 10, T = 16, N = 160 Residuals: Min. 1st Qu. Median 3rd Qu. Max. -154.17714 -11.85856 0.15505 10.77434 142.68974 Coefficients: Estimate Std. Error t-value Pr(>|t|) value 0.073680 0.010987 6.7061 3.945e-10 *** capital 0.204658 0.022589 9.0599 7.158e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Total Sum of Squares: 377500 Residual Sum of Squares: 222510 R-Squared: 0.41058 Adj. R-Squared: 0.36677 F-statistic: 51.5477 on 2 and 148 DF, p-value: < 2.22e-16
変量効果モデルの実行。
こちらも同様にvalueとcapitalの両方が統計的に有意である。
価値valueが1単位増加すると、投資 inv は約0.0813単位増加する
資本capitalが1単位増加すると、投資 inv は約0.2083単位増加する
> #-- 変量効果モデル > fit2 <- plm(inv ~ value + capital, data=ADTS_T, model="random") > summary(fit2) Oneway (individual) effect Random Effect Model (Swamy-Arora's transformation) Call: plm(formula = inv ~ value + capital, data = ADTS_T, model = "random") Balanced Panel: n = 10, T = 16, N = 160 Effects: var std.dev share idiosyncratic 1503.42 38.77 0.212 individual 5596.92 74.81 0.788 theta: 0.8715 Residuals: Min. 1st Qu. Median 3rd Qu. Max. -133.3387 -13.4845 -0.5492 8.9296 158.9204 Coefficients: Estimate Std. Error z-value Pr(>|z|) (Intercept) -5.1914084 26.5668704 -0.1954 0.8451 value 0.0813006 0.0095707 8.4948 <2e-16 *** capital 0.2082550 0.0221473 9.4032 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Total Sum of Squares: 441150 Residual Sum of Squares: 237130 R-Squared: 0.46247 Adj. R-Squared: 0.45562 Chisq: 135.077 on 2 DF, p-value: < 2.22e-16
ハウスマン検定の実行。
今回の場合、p値=0.1931で、帰無仮説が棄却されないため、変量効果モデルが妥当であると判断される。
ただし、ハウスマン検定や変量効果モデルは注意が多い。現実的に相関がない仮定は考えづらく、固定効果モデルを使う方が一般的？
> #-- ハウスマン検定 > phtest(fit1, fit2) Hausman Test data: inv ~ value + capital chisq = 3.289, df = 2, p-value = 0.1931 alternative hypothesis: one model is inconsistent
推定モデルでの予測
1950年までのデータで作成したモデルで、1951年以降を予測してみる。
事前準備として、pdata.frameを使って専用のデータフレームを作成する。index に属性情報と経時情報を指定する。
作成したデータフレームをもとにpredictを使って予測値を算出する。
> #-- 評価指標 > ADTS_P <- pdata.frame(Grunfeld[Grunfeld$FLG==0,], index=c("firm","year")) > pred1 <- predict(fit1, newdata=ADTS_P) > pred2 <- predict(fit2, newdata=ADTS_P)
予測値と実測値でのRMSE、MAE、MAPEを算出。
外れ値の考慮や実測値に0を含むかでどの指標を使うか検討。
> ADTS_P$pred1 <- pred1 > ADTS_P$pred2 <- pred2 > > fnc1 <- function(OBS,PRED){ + RMSE <- sqrt(1/length(PRED) * sum(OBS - PRED)^2) + MAE <- 1/length(PRED) * sum(abs(OBS - PRED)) + MAPE <- 1/length(PRED) * sum(abs(OBS - PRED) / PRED) * 100 + print(paste0("RMSE:",round(RMSE,4)," MAE:",round(MAE,4)," MAPE:",round(MAPE,4))) + } > fnc1(ADTS_P$inv, ADTS_P$pred1) [1] "RMSE:204.1619 MAE:61.3335 MAPE:27.1647" > fnc1(ADTS_P$inv, ADTS_P$pred2) [1] "RMSE:181.0856 MAE:109.2446 MAPE:52.0501"
参考
統計学（出版：東京図書）, 日本統計学会編
https://www.jil.go.jp/institute/zassi/backnumber/2015/04/pdf/006-009.pdf
