#### Hausman test (original version and regression-based version)
##
##
## (1) comparison to Baltagi (2013), sec. 4.3.1, example 1 (pp. 81-82)
## (2) comparison to Baltagi (2013), sec. 4.3.2, example 2 (pp. 82-83)
## (3) comparison to Stata
################################## (1) ##################################
# Baltagi (2013), Econometric Analysis of Panel Data, 5th edition, Wiley & Sons
# Sec 4.3.1, p. 81 (example 1):
#
#### statistics are: 2.33 for original Hausman (m1)
# 2.131 for m2, m3 (for the Grunfeld data)
#
#### vcov within * 10^-3:
#
# 0.14058 -0.077468
# 0.3011788
#
#### vcov between * 10^-3:
#
# 0.82630142 -3.7002477
# 36.4572431
options(digits = 10)
library(plm)
data("Grunfeld", package = "plm")
Grunfeldpdata <- pdata.frame(Grunfeld, index = c("firm", "year"), drop.index = FALSE, row.names = TRUE)
fe_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="within")
be_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="between")
re_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="random")
pool_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="pooling")
# Hausman test
# m1, m2, m3 are all mathematically identical; however computer computation differs a little bit
phtest(inv ~ value + capital, Grunfeldpdata) # replicates Baltagi's m1 = 2.33
phtest(fe_grun, re_grun) # same as above, replicates Baltagi's m1 = 2.33
phtest(re_grun, fe_grun)
phtest(be_grun, re_grun) # replicates Baltagi's m2 = 2.131
phtest(re_grun, be_grun)
phtest(be_grun, fe_grun) # replicates Baltagi's m3 = 2.131 [values m2 and m3 coincide in this case]
phtest(fe_grun, be_grun)
phtest(inv ~ value + capital, Grunfeldpdata, method="aux") # replicates m3 from above in regression test
phtest(inv ~ value + capital, Grunfeldpdata, method="aux", vcov = vcovHC) # no comparison value given
# replicates variance-covariance matrices
vcov(fe_grun)*1000
vcov(be_grun)*1000
################################## (2) ##################################
# Baltagi (2013), Econometric Analysis of Panel Data, 5th edition, Wiley & Sons
# Sec 4.3.2, p. 82-83 (example 2):
### Baltagi's Gasoline example
data("Gasoline", package = "plm")
form <- lgaspcar ~ lincomep + lrpmg + lcarpcap
fe <- plm(form, data = Gasoline, model = "within")
be <- plm(form, data = Gasoline, model = "between")
re <- plm(form, data = Gasoline, model = "random")
phtest(fe, re) # replicates Baltagi's m1 = 302.8
phtest(form, data = Gasoline) # same as above (m1)
phtest(be, re) # replicates Baltagi's m2 = 27.45
phtest(be, fe) # replicates Baltagi's m3 = 26.507 almost
phtest(form, data = Gasoline, method = "aux") # chisq = 26.495054, replicates _almost_ Baltagi's m3 = 26.507
# replicates variance-covariance matrices
#
# vcov in Baltagi within:
# 0.539 0.029 -0.205
# 0.194 0.009
# 0.088
#
# vcov in Baltagi between:
# 2.422 -1.694 -1.056
# 1.766 0.883
# 0.680
vcov(fe)*100
vcov(be)*100
##### twoways case ###
fe2_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="within", effect = "twoways")
# be_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="between")
# RE gives warning due to neg. variance estimation
re2_grun <- plm(inv ~ value + capital, data=Grunfeldpdata, model="random", effect = "twoways")
phtest(fe2_grun, re2_grun) # 13.460, p = 0.00194496 [also given by EViews 9.5;
# Baltagi (2013), p. 85 has other values due to older/wrong version of EViews?]
phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "twoways")
phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "time")
# test to see of phtest(, method = "aux") respects argument effect
# formal test (statistic is about 13 for twoways case and well below in one-way cases)
testobj <- phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "twoways", method = "aux")
testobj2 <- phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "twoways") # just to be sure: test for method="chisq" also...
# test for class of statistic [was matrix pre rev. 305]
testobj1 <- phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "twoways", method = "aux")
testobj2 <- phtest(fe2_grun, re2_grun)
testobj3 <- phtest(inv ~ value + capital, data=Grunfeldpdata, effect = "twoways")
if (class(testobj1$statistic) != "numeric") stop(paste0("class of statistic is not numeric, but ", class(testobj1$statistic)))
if (class(testobj2$statistic) != "numeric") stop(paste0("class of statistic is not numeric, but ", class(testobj2$statistic)))
if (class(testobj3$statistic) != "numeric") stop(paste0("class of statistic is not numeric, but ", class(testobj3$statistic)))
# Two-ways case with between model should result in informative errors.
# phtest(fe2_grun, be_grun)
# phtest(re2_grun, be_grun)
################################## (3) ##################################
### comparison to Stata:
# Hausman test with Stata example 2, pp. 5-6 in http://www.stata.com/manuals/xtxtregpostestimation.pdf
#
# Results of phtest differ, most likely because RE model differs slightly from Stata's RE model as the
# default RE model in Stata uses a slightly different implementation of Swamy-Arora method
# [see http://www.stata.com/manuals/xtxtreg.pdf]
#
# Stata:
# chi2(8) = (b-B)'[(V_b-V_B)^(-1)](b-B)
# = 149.43
# Prob>chi2 = 0.0000
#
# library(haven)
# nlswork <- read_dta("http://www.stata-press.com/data/r14/nlswork.dta") # large file
# nlswork$race <- factor(nlswork$race) # convert
# nlswork$race2 <- factor(ifelse(nlswork$race == 2, 1, 0)) # need this variable for example
# nlswork$grade <- as.numeric(nlswork$grade)
# nlswork$age2 <- (nlswork$age)^2
# nlswork$tenure2 <- (nlswork$tenure)^2
# nlswork$ttl_exp2 <- (nlswork$ttl_exp)^2
#
# pnlswork <- pdata.frame(nlswork, index=c("idcode", "year"), drop.index=F)
#
# form_nls_ex2 <- formula(ln_wage ~ grade + age + age2 + ttl_exp + ttl_exp2 + tenure + tenure2 + race2 + not_smsa + south)
#
# plm_fe_nlswork <- plm(form_nls_ex2, data = pnlswork, model = "within")
# plm_be_nlswork <- plm(form_nls_ex2, data = pnlswork, model = "between")
# plm_re_nlswork <- plm(form_nls_ex2, data = pnlswork, model = "random")
# plm_re2_nlswork <- plm(form_nls_ex2, data = pnlswork, model = "random", random.models = c("within", "between"))
#
# # both slightly different to Stata's example due to diff. RE estimation in Stata
# summary(plm_re_nlswork)
# summary(plm_re2_nlswork)
#
# ### Stata: chi2(8) = 149.43
# phtest(plm_fe_nlswork, plm_re_nlswork) # chisq = 176.39, df = 8, p-value < 2.2e-16
# phtest(plm_be_nlswork, plm_re_nlswork) # chisq = 141.97, df = 10, p-value < 2.2e-16
# phtest(form_nls_ex2, data = pnlswork, method="aux") # chisq = 627.46, df = 8, p-value < 2.2e-16
# phtest(form_nls_ex2, data = nlswork, method="aux") # same on data.frame
# phtest(form_nls_ex2, data = pnlswork, method="aux", vcov = vcovHC) # chisq = 583.56, df = 8, p-value < 2.2e-16
# phtest(form_nls_ex2, data = pnlswork, method="aux", vcov = function(x) vcovHC(x, method="white2", type="HC3"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.