使用するデータ

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