This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
library(plm)
library(lmtest)
library(car)
library(stargazer)data("Grunfeld", package = "plm")pdata <- pdata.frame(Grunfeld, index = c("firm", "year"))pdim(pdata)Balanced Panel: n = 10, T = 20, N = 200summary(pdata) firm year inv value capital
1 :20 1935 : 10 Min. : 0.93 Min. : 58.12 Min. : 0.80
2 :20 1936 : 10 1st Qu.: 33.56 1st Qu.: 199.97 1st Qu.: 79.17
3 :20 1937 : 10 Median : 57.48 Median : 517.95 Median : 205.60
4 :20 1938 : 10 Mean : 145.96 Mean :1081.68 Mean : 276.02
5 :20 1939 : 10 3rd Qu.: 138.04 3rd Qu.:1679.85 3rd Qu.: 358.10
6 :20 1940 : 10 Max. :1486.70 Max. :6241.70 Max. :2226.30
(Other):80 (Other):140 form <- inv ~ value + capitalcem <- plm(form, data = pdata, model = "pooling")fem <- plm(form, data = pdata, model = "within")rem <- plm(form, data = pdata, model = "random")summary(cem)Pooling Model
Call:
plm(formula = form, data = pdata, model = "pooling")
Balanced Panel: n = 10, T = 20, N = 200
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-291.6757 -30.0137 5.3033 34.8293 369.4464
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) -42.7143694 9.5116760 -4.4907 1.207e-05 ***
value 0.1155622 0.0058357 19.8026 < 2.2e-16 ***
capital 0.2306785 0.0254758 9.0548 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 9359900
Residual Sum of Squares: 1755900
R-Squared: 0.81241
Adj. R-Squared: 0.8105
F-statistic: 426.576 on 2 and 197 DF, p-value: < 2.22e-16summary(fem)Oneway (individual) effect Within Model
Call:
plm(formula = form, data = pdata, model = "within")
Balanced Panel: n = 10, T = 20, N = 200
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-184.00857 -17.64316 0.56337 19.19222 250.70974
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
value 0.110124 0.011857 9.2879 < 2.2e-16 ***
capital 0.310065 0.017355 17.8666 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 2244400
Residual Sum of Squares: 523480
R-Squared: 0.76676
Adj. R-Squared: 0.75311
F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16summary(rem)Oneway (individual) effect Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = form, data = pdata, model = "random")
Balanced Panel: n = 10, T = 20, N = 200
Effects:
var std.dev share
idiosyncratic 2784.46 52.77 0.282
individual 7089.80 84.20 0.718
theta: 0.8612
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-177.6063 -19.7350 4.6851 19.5105 252.8743
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) -57.834415 28.898935 -2.0013 0.04536 *
value 0.109781 0.010493 10.4627 < 2e-16 ***
capital 0.308113 0.017180 17.9339 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 2381400
Residual Sum of Squares: 548900
R-Squared: 0.7695
Adj. R-Squared: 0.76716
Chisq: 657.674 on 2 DF, p-value: < 2.22e-16fixef(fem) 1 2 3 4 5 6 7 8 9 10
-70.2967 101.9058 -235.5718 -27.8093 -114.6168 -23.1613 -66.5535 -57.5457 -87.2223 -6.5678 stargazer(cem, fem, rem, type = "text", column.labels = c("CEM", "FEM", "REM"))Warning: length of NULL cannot be changedWarning: length of NULL cannot be changedWarning: length of NULL cannot be changedWarning: length of NULL cannot be changedWarning: length of NULL cannot be changedWarning: number of rows of result is not a multiple of vector length (arg 2)
=========================================================================
Dependent variable:
------------------------------------------------------------
inv
CEM FEM REM
(1) (2) (3)
-------------------------------------------------------------------------
value 0.116*** 0.110*** 0.110***
(0.006) (0.012) (0.010)
capital 0.231*** 0.310*** 0.308***
(0.025) (0.017) (0.017)
Constant -42.714*** -57.834**
(9.512) (28.899)
-------------------------------------------------------------------------
Observations 200 200 200
R2 0.812 0.767 0.770
Adjusted R2 0.811 0.753 0.767
F Statistic 426.576*** (df = 2; 197) 309.014*** (df = 2; 188) 657.674***
=========================================================================
Note: *p<0.1; **p<0.05; ***p<0.01pFtest(fem, cem)
F test for individual effects
data: form
F = 49.177, df1 = 9, df2 = 188, p-value < 2.2e-16
alternative hypothesis: significant effectsplmtest(cem, type = "bp")
Lagrange Multiplier Test - (Breusch-Pagan)
data: form
chisq = 798.16, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effectsphtest(fem, rem)
Hausman Test
data: form
chisq = 2.3304, df = 2, p-value = 0.3119
alternative hypothesis: one model is inconsistentfem_tw <- plm(form, data = pdata, model = "within", effect = "twoways")
pFtest(fem_tw, fem)
F test for twoways effects
data: form
F = 1.4032, df1 = 19, df2 = 169, p-value = 0.1309
alternative hypothesis: significant effectsbptest(inv ~ value + capital + factor(firm), data = Grunfeld)
studentized Breusch-Pagan test
data: inv ~ value + capital + factor(firm)
BP = 85.836, df = 11, p-value = 1.086e-13pbgtest(fem)
Breusch-Godfrey/Wooldridge test for serial correlation in panel models
data: form
chisq = 65.063, df = 20, p-value = 1.14e-06
alternative hypothesis: serial correlation in idiosyncratic errorspwartest(fem)
Wooldridge's test for serial correlation in FE panels
data: fem
F = 76.929, df1 = 1, df2 = 188, p-value = 1.057e-15
alternative hypothesis: serial correlationpcdtest(fem, test = "cd")
Pesaran CD test for cross-sectional dependence in panels
data: inv ~ value + capital
z = 4.6612, p-value = 3.144e-06
alternative hypothesis: cross-sectional dependencevif(lm(inv ~ value + capital, data = Grunfeld)) value capital
1.313773 1.313773 cor(Grunfeld[, c("value", "capital")]) value capital
value 1.0000000 0.4887058
capital 0.4887058 1.0000000purtest(pdata$inv, test = "levinlin", exo = "intercept", lags = "AIC", pmax = 4)Warning: the time series is short
Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
data: pdata$inv
z = 3.1416, p-value = 0.9992
alternative hypothesis: stationaritycoeftest(fem, vcov = vcovHC(fem, type = "HC1", cluster = "group"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
value 0.110124 0.014414 7.6398 1.076e-12 ***
capital 0.310065 0.050043 6.1959 3.576e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1coeftest(fem, vcov = vcovSCC(fem, type = "HC1"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
value 0.110124 0.017775 6.1954 3.586e-09 ***
capital 0.310065 0.034996 8.8601 6.017e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1