1 2R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out" 3Copyright (C) 2020 The R Foundation for Statistical Computing 4Platform: x86_64-w64-mingw32/x64 (64-bit) 5 6R is free software and comes with ABSOLUTELY NO WARRANTY. 7You are welcome to redistribute it under certain conditions. 8Type 'license()' or 'licence()' for distribution details. 9 10R is a collaborative project with many contributors. 11Type 'contributors()' for more information and 12'citation()' on how to cite R or R packages in publications. 13 14Type 'demo()' for some demos, 'help()' for on-line help, or 15'help.start()' for an HTML browser interface to help. 16Type 'q()' to quit R. 17 18> # test pbsytest() - unbalanced and balanced version 19> 20> ################### Bera, Sosa-Escudero and Yoon (2001) and joint test of Baltagi/Li (1991) ############### 21> # see Baltagi (2005), Econometric Analysis of Panel Data, 3rd edition, pp. 96-97 22> # or Baltagi (2013), Econometric Analysis of Panel Data, 5th edition, p. 108. 23> # 24> ## only balanced tests described in Bera, Sosa-Escudero and Yoon (2001) and Baltagi (2005, 2013)! 25> # 26> # Baltagi (2013), p. 108: 27> # Grunfeld data, (table 4.2) 28> # LM_mu = 798.162 (with Stata's xttest0 command) [-> plmtest(pool_grunfeld, type = "bp")] 29> # LM_rho = 143.523, LM*_mu = 664.948, LM*_rho = 10.310, joint test (LM1) = 808.471 (all using TSP) 30> # 31> # comments about significance in book: 32> # joint test (LM1): rejects null hypo (no first-order serial correlation and no random effects) 33> # LM_rho, LM*_rho: reject null hypo (no first-order serial correlation) 34> # LM_mu, LM*_mu: reject null hypo (no random effects) 35> 36> 37> library(plm) 38> data("Grunfeld", package = "plm") 39> Grunfeldpdata <- pdata.frame(Grunfeld, index = c("firm", "year"), drop.index = FALSE, row.names = TRUE) 40> form_gunfeld <- formula(inv ~ value + capital) 41> pool_grunfeld <- plm(form_gunfeld, data = Grunfeldpdata, model="pooling") 42> 43> pbsytest(pool_grunfeld, test = "ar") # chisq = 10.31 => LM*_rho in Baltagi's book (RS*_lambda from Sosa-Escudero/Bera (2008), p. 73) 44 45 Bera, Sosa-Escudero and Yoon locally robust test - balanced panel 46 47data: formula 48chisq = 10.31, df = 1, p-value = 0.001323 49alternative hypothesis: AR(1) errors sub random effects 50 51> pbsytest(pool_grunfeld, test = "re", re.normal = FALSE) # chisq = 664.948 => LM*_mu in Baltagi's book (RS*_mu from Sosa-Escudero/Bera (2008), p. 73) 52 53 Bera, Sosa-Escudero and Yoon locally robust test (two-sided) - 54 balanced panel 55 56data: formula 57chisq = 664.95, df = 1, p-value < 2.2e-16 58alternative hypothesis: random effects sub AR(1) errors 59 60> pbsytest(pool_grunfeld, test = "re") # [sqrt(chisq) = z = 25.787] => RSO*_mu from Sosa-Escudero/Bera (2008), p. 75 61 62 Bera, Sosa-Escudero and Yoon locally robust test (one-sided) - 63 balanced panel 64 65data: formula 66z = 25.787, p-value < 2.2e-16 67alternative hypothesis: random effects sub AR(1) errors 68 69> pbsytest(pool_grunfeld, test = "j") # chisq = 808.47 => LM1 in Baltagi's book (RS_lambda_mu in Sosa-Escudero/Bera (2008), p. 74) 70 71 Baltagi and Li AR-RE joint test - balanced panel 72 73data: formula 74chisq = 808.47, df = 2, p-value < 2.2e-16 75alternative hypothesis: AR(1) errors or random effects 76 77> 78> # formula interface 79> pbsytest(form_gunfeld, data = Grunfeld, test = "ar") 80 81 Bera, Sosa-Escudero and Yoon locally robust test - balanced panel 82 83data: formula 84chisq = 10.31, df = 1, p-value = 0.001323 85alternative hypothesis: AR(1) errors sub random effects 86 87> pbsytest(form_gunfeld, data = Grunfeld, test = "re") 88 89 Bera, Sosa-Escudero and Yoon locally robust test (one-sided) - 90 balanced panel 91 92data: formula 93z = 25.787, p-value < 2.2e-16 94alternative hypothesis: random effects sub AR(1) errors 95 96> pbsytest(form_gunfeld, data = Grunfeld, test = "re", re.normal = FALSE) 97 98 Bera, Sosa-Escudero and Yoon locally robust test (two-sided) - 99 balanced panel 100 101data: formula 102chisq = 664.95, df = 1, p-value < 2.2e-16 103alternative hypothesis: random effects sub AR(1) errors 104 105> pbsytest(form_gunfeld, data = Grunfeld, test = "j") 106 107 Baltagi and Li AR-RE joint test - balanced panel 108 109data: formula 110chisq = 808.47, df = 2, p-value < 2.2e-16 111alternative hypothesis: AR(1) errors or random effects 112 113> 114> plmtest(pool_grunfeld, type = "bp") # LM_mu in Baltagi's book 115 116 Lagrange Multiplier Test - (Breusch-Pagan) for balanced panels 117 118data: form_gunfeld 119chisq = 798.16, df = 1, p-value < 2.2e-16 120alternative hypothesis: significant effects 121 122> 123> ############### balanced version ################### 124> ### Results from Bera et al. (2001), p. 13: 125> 126> ## Bera/Sosa-Escudero/Yoon (2001), Tests for the error component model in the presence of local misspecifcation, 127> ## Journal of Econometrics 101 (2001), pp. 1-23. 128> 129> # To replicate, a special version of the Grunfeld data set is needed: only 5 selected firms (total of 100 obs) 130> # from http://pages.stern.nyu.edu/~wgreene/Text/tables/TableF13-1.txt 131> # or http://statmath.wu.ac.at/~zeileis/grunfeld/TableF13-1.txt 132> # 133> # NB: this data set contains 3 errors compared to the original Grunfeld data, see e.g., the 134> # analysis of various different Grundfeld data sets circulating at http://statmath.wu-wien.ac.at/~zeileis/grunfeld/ 135> # or https://eeecon.uibk.ac.at/~zeileis/grunfeld/ 136> # 137> ## commented due to file download 138> 139> # Grunfeld_greene_5firms <- read.csv("http://pages.stern.nyu.edu/~wgreene/Text/tables/TableF13-1.txt", sep="") 140> # # Grunfeld_greene_5firms <- read.csv("http://statmath.wu.ac.at/~zeileis/grunfeld/TableF13-1.txt", sep="") # alternative source 141> # 142> # # Matching to Grunfeld data set in plm 143> # # Grunfeld[c(1:20, 41:60), 3:5] == Grunfeld_greene_5firms[c(1:20, 41:60), 3:5] 144> # # Grunfeld[61:80, 3:5] == Grunfeld_greene_5firms[21:40, 3:5] 145> # # Grunfeld[141:160, 3:5] == Grunfeld_greene_5firms[61:80, 3:5] 146> # # Grunfeld[21:40, 3:5] == Grunfeld_greene_5firms[81:100, 3:5] # almost all equal, 3 values differ (3 errors in the Greene 5 firm version) 147> # 148> # pGrunfeld_greene_5firms <- pdata.frame(Grunfeld_greene_5firms, index = c("Firm", "Year"), drop.index = FALSE, row.names = TRUE) 149> # form_gunfeld_half <- formula(I ~ F + C) 150> # pool_grunfeld_half <- plm(form_gunfeld_half, data=pGrunfeld_greene_5firms, model = "pooling") 151> # re_grunfeld_half <- plm(form_gunfeld_half, data=pGrunfeld_greene_5firms, model = "random") 152> # 153> # pbsytest(pool_grunfeld_half, test = "ar") # chisq = 3.7125 => RS*_rho in Bera et al. (2001), p. 13 154> # pbsytest(pool_grunfeld_half, test = "re") # normal = 19.601; p = 0 => RSO*_mu 155> # pbsytest(pool_grunfeld_half, test = "re", re.normal = FALSE) # chisq = 384.183 => RS*_mu [sqrt(chisq) = z = 19.601] 156> # pbsytest(pool_grunfeld_half, test = "j") # chisq = 457.53 => RS_mu_rho 157> # 158> # # plmtest's statistic is also mentioned in paper 159> # plmtest(pool_grunfeld_half, type = "bp") # chisq = 453.82 => RS_mu in Bera et al. (2001), p. 13 160> # plmtest(pool_grunfeld_half, type = "honda") # normal = 21.3031 => RSO_mu 161> # 162> # 163> # ## RS_rho in Bera et al (2001), p. 9 (formula 19) is not implemented 164> # ## it's origin is in Baltagi/Li (1991), but there is is just a side result 165> # ## in terms of n, t, b of pbsystest it is: (n*t^2*(B^2)) / (t-1) 166> # 167> # # formula interface 168> # pbsytest(form_gunfeld_half, data = pGrunfeld_greene_5firms, test = "ar") 169> # pbsytest(form_gunfeld_half, data = pGrunfeld_greene_5firms, test = "re") 170> # pbsytest(form_gunfeld_half, data = pGrunfeld_greene_5firms, test = "re", re.normal = FALSE) 171> # pbsytest(form_gunfeld_half, data = pGrunfeld_greene_5firms, test = "j") 172> # 173> # plmtest(form_gunfeld_half, data = pGrunfeld_greene_5firms, type = "bp") 174> 175> 176> ############ Replicate tests from original paper Sosa-Escudero/Bera (2008) #################### 177> ############ unbalanced panel #################### 178> ## 179> ## data set for test from Sosa-Escudero/Bera (2008), pp. 75-77 180> ## available as Stata .dta file at http://www.stata-journal.com/software/sj8-1/sg164_1/ginipanel5.dta 181> ## 182> ## Sosa-Escudero/Bera (2008), Tests for unbalanced error-components models under local misspecification, 183> ## The Stata Journal (2008), Vol. 8, Number 1, pp. 68-78. 184> 185> ## Commented due to extra package needed 186> 187> # library(haven) 188> # ginipanel5 <- read_dta("http://www.stata-journal.com/software/sj8-1/sg164_1/ginipanel5.dta") 189> # pginipanel5 <- pdata.frame(ginipanel5, index = c("naglo", "ano"), drop.index = FALSE, row.names = TRUE) 190> # 191> # # Stata command for RE model: xtreg gini ie ie2 indus adpubedsal desempleo tactiv invipib apertura pyas4 e64 supc tamfam, re 192> # # use pooling model in R: 193> # formula_gini <- formula(gini ~ ie + ie2 + indus + adpubedsal + desempleo + tactiv + invipib + apertura + pyas4 + e64 + supc + tamfam) 194> # pool_gini <- plm(formula_gini, data = pginipanel5, model = "pooling") 195> # 196> # pdim(pool_gini) # Unbalanced Panel: n=17, T=6-8, N=128 197> # 198> # # Stata's Output of xttest1, unadjusted (Sosa-Escudero/Bera (2008), p. 77): 199> # # 200> # # Random Effects, Two Sided: 201> # # LM(Var(u)=0) = 13.50 Pr>chi2(1) = 0.0002 202> # # ALM(Var(u)=0) = 6.03 Pr>chi2(1) = 0.0141 # test="re", re.normal = FALSE 203> # # 204> # # Random Effects, One Sided: 205> # # LM(Var(u)=0) = 3.67 Pr>N(0,1) = 0.0001 206> # # ALM(Var(u)=0) = 2.46 Pr>N(0,1) = 0.0070 # test="re", re.normal = TRUE 207> # # 208> # # Serial Correlation: 209> # # LM(lambda=0) = 9.32 Pr>chi2(1) = 0.0023 210> # # ALM(lambda=0) = 1.86 Pr>chi2(1) = 0.1732 # test="ar" 211> # # 212> # # Joint Test: 213> # # LM(Var(u)=0,lambda=0) = 15.35 Pr>chi2(2) = 0.0005 # test="j" 214> # 215> # 216> # pbsytest(pool_gini, test = "re", re.normal = FALSE) # chisq = 6.0288793, df = 1, p-value = 0.01407367 217> # pbsytest(pool_gini, test = "re") # normal = 2.4553776, n/a p-value = 0.007036833 218> # pbsytest(pool_gini, test = "ar") # chisq = 1.8550073, df = 1, p-value = 0.1732021 219> # pbsytest(pool_gini, test = "j") # chisq = 15.352307, df = 2, p-value = 0.0004637552 220> # 221> # # formula interface 222> # pbsytest(formula_gini, data = pginipanel5, test = "re", re.normal = FALSE) # chisq = 6.0288793, df = 1, p-value = 0.01407367 223> # pbsytest(formula_gini, data = pginipanel5, test = "re") # normal = 2.4553776, n/a p-value = 0.007036833 224> # pbsytest(formula_gini, data = pginipanel5, test = "ar") # chisq = 1.8550073, df = 1, p-value = 0.1732021 225> # pbsytest(formula_gini, data = pginipanel5, test = "j") # chisq = 15.352307, df = 2, p-value = 0.0004637552 226> # 227> 228> proc.time() 229 user system elapsed 230 0.62 0.17 0.78 231