Changes between Initial Version and Version 1 of PanelLinearModel01

Timestamp:

Nov 14, 2016, 10:15:31 AM (7 years ago)

Author:

admin

Comment:

--

PanelLinearModel01

@@ +1 @@
+R言語による多変量時系列分析。
+
+複数グループ・複数項目の時系列データで時系列間の関係性・影響を明らかにする。
+
+
+以下、実行可能なR言語ソースコードも用い紹介する。
+
+
+例：複数時系列間の関係性・影響
+
+農家ごとの農地の肥沃度・人の各仕事の労働量・各肥料の投入量の各時系列、および、生産量の時系列から、肥沃度・各労働量・各肥料量が生産量に与える影響を明らかにする。
+国ごとの各政策種類への投資額の時系列、および、GDP から、各政策投資がGDPに与える影響を明らかにする
+地域ごとの各キャンペーン種類の投入量の各時系列、および、商品売上の時系列から、各地域の各キャンペーンが売上に与える影響を明らかにする。等。
+
+R Library
+
+複数時系列の時系列分析として、R言語では plm 、 ccgarch 等のパッケージがあるが、今回は plm (Panel Linear Model) を紹介する。
+
+※参考：Panel Data Econometrics in R: The plm Package (PDF)
+
+
+Panel Data(交差時系列)
+
+計量経済学の領域では、交差系列 (Cross section Data) で時系列 (Time series Data) である系列を
+
+『交差時系列(Panel Data)』と呼ぶ。　
+
+{{{
+例）47都道府県の人口を時間を追って調べたもの。
+
+交差系列 (Cross section Data)：同一時点での様々なData。例) ある時点で47都道府県の人口、人口密度、男女比などを調べたもの。
+時系列 (Time series Data)：同一種類のDataを様々な時点で取ったもの。例) ある都道府県の人口を時間を追って調べたもの。
+※参考：Wikipedia: 計量経済学 - Wikipedia、パネルデータの意義とその活用(PDF)
+
+
+関数
+
+plm package には複数時系列の評価関数として次の4つが用意されている。
+
+
+Functions:
+
+plm: estimation of the basic panel models, i.e. within, between and random e
+ffect models. Models are estimated using the lm function to transformed data,
+pvcm: estimation of models with variable coecients,
+pgmm: estimation of generalized method of moments models,
+pggls: estimation of general feasible generalized least squares models
+これらはRの線形回帰関数 lm() と同じ引数の形をとり、1引数目に formula、2引数目にdata をとる。また、Option として、次の３つを指定できる。
+
+
+Options:
+
+index: this argument enables the estimation functions to identify the structure of the data, i.e. the individual and the time period for each observation,
+effect: the kind of e
+ffects to include in the model, i.e. individual e
+ffects, time e
+ffects or both,
+model: the kind of model to be estimated, most of the time a model with 
+xed e
+ffects or a model with random e
+ffects.
+
+Results:
+
+上記各関数の結果は 関数名と同じクラスのオブジェクトに保存され、これらのクラスは panelmodel オブジェクトを継承している。これらのオブジェクトは以下を含んでいる。
+
+coefficients
+residuals
+fitted.values
+vcov
+df.residual
+call
+
+データ
+
+Rの組み込みデータ「Produc： 米国の48州の1970-1986のパネルデータ」を用いる
+
+library("plm")
+data("Produc", package = "plm")
+head(Produc, 10)
+Produc データ (10行目まで抜粋)：
+
+> head(Produc, 10)
+     state year     pcap     hwy   water    util       pc   gsp    emp unemp
+1  ALABAMA 1970 15032.67 7325.80 1655.68 6051.20 35793.80 28418 1010.5   4.7
+2  ALABAMA 1971 15501.94 7525.94 1721.02 6254.98 37299.91 29375 1021.9   5.2
+3  ALABAMA 1972 15972.41 7765.42 1764.75 6442.23 38670.30 31303 1072.3   4.7
+4  ALABAMA 1973 16406.26 7907.66 1742.41 6756.19 40084.01 33430 1135.5   3.9
+5  ALABAMA 1974 16762.67 8025.52 1734.85 7002.29 42057.31 33749 1169.8   5.5
+6  ALABAMA 1975 17316.26 8158.23 1752.27 7405.76 43971.71 33604 1155.4   7.7
+7  ALABAMA 1976 17732.86 8228.19 1799.74 7704.93 50221.57 35764 1207.0   6.8
+8  ALABAMA 1977 18111.93 8365.67 1845.11 7901.15 51084.99 37463 1269.2   7.4
+9  ALABAMA 1978 18479.74 8510.64 1960.51 8008.59 52604.05 39964 1336.5   6.3
+10 ALABAMA 1979 18881.49 8640.61 2081.91 8158.97 54525.86 40979 1362.0   7.1
+state（州）、year(年）、gsp (州内総生産）、pcap (民間資本ストック）、pc(公的資本）、emp (労働投入量）
+
+※参考：Rパッケージガイドブック 蓮見亮さん記事。
+
+
+Rでの実装例
+
+library("plm")
+
+# Producデータ
+data("Produc", package = "plm")
+data <- Produc
+
+# index に用いる グループ列名、時間列名
+grpCName <- "state"
+timeCName <- "year"
+
+# pdata.frame (panel data frame) 作成。グループ列, 時間列で index 作成　 
+idx <- c(grpCName, timeCName)
+pdata <- pdata.frame(data, index = idx, drop.index = TRUE, row.names = TRUE)
+
+# 目的変数列 指定
+tidx <- 6
+
+###################################
+# Panel Linear Model 作成
+###################################
+# plm(formula, data, subset, na.action, effect = c("individual","time","twoways"),
+#    model = c("within","random","ht","between","pooling","fd"),
+#    random.method = c("swar","walhus","amemiya","nerlove"),
+#    inst.method = c("bvk","baltagi"), index = NULL, ...)
+#
+#[model 指定]：
+# within: Fixed Effect Model [Default] 
+#   (固定効果モデル。個別効果を完全除去。
+#    目的変数、説明変数、残差のそれぞれをグループごとの平均を回帰式から引き算)
+# pooling: Pooling Model     (グループ、時間を気にせず回帰）
+# fd: First-Difference Model (グループごとに回帰。変数の時間差分を用いる）
+# between: Between Model　　　　（グループごとに回帰。変数の時間平均を用いる）
+# random: error components model
+#   (変量効果モデル。残差をグループ間のものとグループ内のものに分割）。
+#    ramdom.method 変動性の指定ができ swar [Default], walhus, amemiya, nerlove を指定可能。
+#[effect 指定]:
+# effect = "time" : 上記処理をグループ毎の処理に置き換える (グループと時間を入れ替え）
+
+# Formula作成：目的変数以外を説明変数とする
+cn <- colnames(pdata)
+fmla <- as.formula(paste(cn[tidx], " ~ ", paste(cn[-tidx], collapse = " + ")))
+
+# 固定効果モデル
+res.fe <-  plm(formula = fmla, data = pdata, model = "within")
+#(res.fe <-  plm(formula = fmla, data = data, index = idx, model = "within"))　と同じ
+summary(res.fe)
+
+# 変動効果モデル
+res.re <-  plm(formula = fmla, data = pdata, model = "random")
+summary(res.re)
+
+###################################
+# 各種検定
+###################################
+
+# 各種豊富な検定関数が用意されている。
+
+library(car)
+# Tests for individual and time effect
+plmtest(fmla, data = pdata, effect = "twoways", type = "ghm")
+# Tests of poolability
+pooltest(fmla, data = pdata, model = "within")
+# F tests of e^Kffects based on the comparison of the within and the pooling models
+pFtest(fmla, data = pdata, effect = "twoways")
+# Hausman test based on the comparison of two sets of estimates
+phtest(res.fe, res.re)
+# Tests of serial correlatio
+pwtest(fmla, data = pdata)
+# Locally robust tests for serial correlation or random effects
+pbsytest(fmla, data = pdata, test = "j")
+#Conditional LM test for AR(1) or MA(1) errors under random effect
+pbltest(fmla, data = data, alternative = "onesided") #pdata.frame ではなく元 data.frame を用いる
+
+library(lmtest)
+#General serial correlation test
+pbgtest(res.fe, order = 2)
+#Wooldridge's test for serial correlation in short FE panels
+pwartest(fmla, data = pdata)
+#Wooldridge's first-difference-based test
+pwfdtest(fmla, data = pdata)
+
+# 他にも Cross-sectional dependence の検定も用意されている。
+
+実行結果
+
+残差、回帰係数、R^2、等、各種回帰結果が得られている。
+
+> # 固定効果モデル
+> res.fe <-  plm(formula = fmla, data = pdata, model = "within")
+> summary(res.fe)
+Oneway (individual) effect Within Model
+
+Call:
+plm(formula = fmla, data = pdata, model = "within")
+
+Balanced Panel: n=48, T=17, N=816
+
+Residuals :
+     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
+-1.54e+04 -7.66e+02 -7.04e+01  1.64e-09  7.71e+02  1.68e+04 
+
+Coefficients :
+         Estimate  Std. Error t-value  Pr(>|t|)    
+pcap   1.1832e+04  1.5029e+04  0.7873   0.43136    
+hwy   -1.1833e+04  1.5029e+04 -0.7873   0.43134    
+water -1.1831e+04  1.5029e+04 -0.7872   0.43141    
+util  -1.1832e+04  1.5029e+04 -0.7873   0.43134    
+pc     1.2001e-01  2.1647e-02  5.5437 4.086e-08 ***
+emp    3.4932e+01  8.8757e-01 39.3571 < 2.2e-16 ***
+unemp -1.6813e+02  6.0568e+01 -2.7759   0.00564 ** 
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
+
+Total Sum of Squares:    1.3289e+11
+Residual Sum of Squares: 4.978e+09
+R-Squared      :  0.96254 
+      Adj. R-Squared :  0.89766 
+F-statistic: 2793.44 on 7 and 761 DF, p-value: < 2.22e-16
+
+> # 変動効果モデル
+> res.re <-  plm(formula = fmla, data = pdata, model = "random")
+summary(res.re)
+> Oneway (individual) effect Random Effect Model 
+   (Swamy-Arora's transformation)
+
+Call:
+plm(formula = fmla, data = pdata, model = "random")
+
+Balanced Panel: n=48, T=17, N=816
+
+Effects:
+                   var  std.dev share
+idiosyncratic  6541433     2558 0.185
+individual    28839985     5370 0.815
+theta:  0.8853  
+
+Residuals :
+    Min.  1st Qu.   Median  3rd Qu.     Max. 
+-11500.0   -847.0     61.7    849.0  19800.0 
+
+Coefficients :
+               Estimate  Std. Error t-value  Pr(>|t|)    
+(Intercept) -1.7738e+03  1.1556e+03 -1.5350    0.1252    
+pcap         1.4200e+04  1.5615e+04  0.9094    0.3634    
+hwy         -1.4201e+04  1.5615e+04 -0.9094    0.3634    
+water       -1.4199e+04  1.5615e+04 -0.9093    0.3634    
+util        -1.4201e+04  1.5615e+04 -0.9094    0.3634    
+pc           1.3118e-01  1.7636e-02  7.4380 2.613e-13 ***
+emp          3.4150e+01  7.6819e-01 44.4552 < 2.2e-16 ***
+unemp       -2.7338e+02  5.9394e+01 -4.6028 4.840e-06 ***
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
+
+Total Sum of Squares:    1.8368e+11
+Residual Sum of Squares: 5743500000
+R-Squared      :  0.96873 
+      Adj. R-Squared :  0.95923 
+F-statistic: 3576.08 on 7 and 808 DF, p-value: < 2.22e-16
+
+検定結果例
+
+> # Tests for individual and time effect
+> plmtest(fmla, data = pdata, effect = "twoways", type = "ghm")
+
+        Lagrange Multiplier Test - two-ways effects (Gourieroux, Holly and
+        Monfort)
+
+data:  fmla 
+chisq = 3156.201, df = 2, p-value < 2.2e-16
+alternative hypothesis: significant effects 
+
+参考資料：時系列分析入門資料
+
+Tokyo.R で講師をした時系列分析の入門資料です。[Ｒ][データマイニング] R言語による時系列分析
+
+}}}