1library(plm) 2data("Grunfeld", package = "plm") 3Grunfeld_unbal <- Grunfeld[1:199, ] 4# ercomp(plm(inv ~ value, Grunfeld, model = "random")) 5# ercomp(plm(inv ~ value, Grunfeld, model = "random", random.method = "amemiya")) 6# ercomp(plm(inv ~ value + capital, Grunfeld_unbal, model = "random")) 7 8 9# these resulted in errors pre rev. 523 due to missing drop = FALSE 10plm(inv ~ value, Grunfeld_unbal, model = "random", random.method = "amemiya") 11plm(inv ~ value, Grunfeld_unbal, model = "random", random.method = "amemiya", effect = "time") 12 13 14# test case for illegal pseries in pmerge's return value: 15# up to rev. 675, pmerge produced a data.frame with a column declared to be a pseries but with lacking index, 16# and there should be no 'pseries' in the resulting data.frame in first place 17pGrunfeld <- pdata.frame(Grunfeld) 18df_after_pmerge <- plm:::pmerge(pGrunfeld$inv, pGrunfeld$value) 19if (inherits(df_after_pmerge$ind, "pseries") && is.null(attr(df_after_pmerge$ind, "index"))) stop("illegal pseries (no index) produced by pmerge") 20if ("pseries" %in% unlist(lapply(df_after_pmerge, class))) stop("pmerge returned a column with pseries") 21if (!"data.frame" == class(df_after_pmerge)) stop("pmerge did not return a pure data.frame according to class()") 22 23 24# pmodel.response: test case for illegal pseries 25form <- formula(inv ~ value + capital) 26if (!is.pseries(pmodel.response(form, data = pGrunfeld, model = "pooling"))) stop("pmodel.response's return value is not a valid pseries") 27if (!is.pseries(pmodel.response(form, data = pGrunfeld, model = "within"))) stop("pmodel.response's return value is not a valid pseries") 28if (!is.pseries(pmodel.response(form, data = pGrunfeld, model = "Between"))) stop("pmodel.response's return value is not a valid pseries") 29if (!is.pseries(pmodel.response(plm(form, data = pGrunfeld, model = "random")))) stop("pmodel.response's return value is not a valid pseries") 30# for FD and between models, it should be a numeric as a pseries does not make sense due to the data compression 31if (inherits(pmodel.response(form, data = pGrunfeld, model = "fd"), "pseries")) stop("pmodel.response's return value shall not be a pseries for fd models") 32if (inherits(pmodel.response(form, data = pGrunfeld, model = "between"), "pseries")) stop("pmodel.response's return value shall not be a pseries for between models") 33if (plm:::has.index(pmodel.response(plm(form, data = pGrunfeld, model = "fd")))) stop("pmodel.response's return value shall not have an index for fd models") 34if (plm:::has.index(pmodel.response(plm(form, data = pGrunfeld, model = "between")))) stop("pmodel.response's return value shall not have an index for between models") 35 36 37# residuals.plm: test case for illegal pseries 38if (!is.pseries(residuals(plm(form, data = pGrunfeld, model = "pooling")))) stop("residuals.plm's return value is not a valid pseries") 39if (!is.pseries(residuals(plm(form, data = pGrunfeld, model = "within")))) stop("residuals.plm's return value is not a valid pseries") 40if (!is.pseries(residuals(plm(form, data = pGrunfeld, model = "random")))) stop("residuals.plm's return value is not a valid pseries") 41# for FD and between models, it should be a numeric as a pseries does not make sense due to the data compression 42if (inherits(residuals(plm(form, data = pGrunfeld, model = "fd")), "pseries")) stop("residuals.plm's return value shall not be a pseries for fd models") 43if (inherits(residuals(plm(form, data = pGrunfeld, model = "between")), "pseries")) stop("residuals.plm's return value shall not be a pseries for between models") 44if (plm:::has.index(residuals(plm(form, data = pGrunfeld, model = "fd")))) stop("residuals.plm's return value shall not have an index for fd models") 45if (plm:::has.index(residuals(plm(form, data = pGrunfeld, model = "between")))) stop("residuals.plm's return value shall not have an index for between models") 46 47 48# fitted.plm: test case for illegal pseries 49if (!is.pseries(fitted(plm(form, data = pGrunfeld, model = "pooling")))) stop("fitted.plm's return value is not a valid pseries") 50if (!is.pseries(fitted(plm(form, data = pGrunfeld, model = "within")))) stop("fitted.plm's return value is not a valid pseries") 51if (!is.pseries(fitted(plm(form, data = pGrunfeld, model = "random")))) stop("fitted.plm's return value is not a valid pseries") 52# for FD and between models, it should be a numeric as a pseries does not make sense due to the data compression 53if (inherits(fitted(plm(form, data = pGrunfeld, model = "fd")), "pseries")) stop("fitted.plm's return value shall not be a pseries for fd models") 54if (inherits(fitted(plm(form, data = pGrunfeld, model = "between")), "pseries")) stop("fitted.plm's return value shall not be a pseries for between models") 55if (plm:::has.index(fitted(plm(form, data = pGrunfeld, model = "fd")))) stop("fitted.plm's return value shall not have an index for fd models") 56if (plm:::has.index(fitted(plm(form, data = pGrunfeld, model = "between")))) stop("fitted.plm's return value shall not have an index for between models") 57 58## WLS 59p <- plm(inv ~ value, Grunfeld, model = "pooling") 60pwls <- plm(inv ~ value + capital, data = Grunfeld, weights = Grunfeld$capital, model = "pooling") 61 62if (!is.null(p$weights)) stop("element 'weights' in plm object albeit it should not be there") 63if (is.null(pwls$weights)) stop("element 'weights' missing in plm object") 64 65## aneweytest 66data("RiceFarms", package = "plm") 67aneweytest(log(goutput) ~ log(seed) + log(totlabor) + log(size), RiceFarms, index = "id") 68 69## piest 70pirice <- piest(log(goutput) ~ log(seed) + log(totlabor) + log(size), RiceFarms, index = "id") 71summary(pirice) 72 73## mtest 74data("EmplUK", package = "plm") 75ar <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) + 76 lag(log(capital), 0:2) + lag(log(output), 0:2) | lag(log(emp), 2:99), 77 data = EmplUK, effect = "twoways", model = "twosteps") 78mtest(ar, order = 1) 79mtest(ar, order = 2, vcov = vcovHC) 80 81## pcdtest 82pcdtest(inv ~ value + capital, data = Grunfeld, 83 index = c("firm", "year")) 84 85## test on two-way fixed effects homogeneous model 86pcdtest(inv ~ value + capital, data = Grunfeld, model = "within", 87 effect = "twoways", index = c("firm", "year")) 88 89## test on panelmodel object 90g <- plm(inv ~ value + capital, data = Grunfeld, index = c("firm", "year")) 91pcdtest(g) 92 93## scaled LM test 94pcdtest(g, test = "sclm") 95 96## test on pseries 97pGrunfeld <- pdata.frame(Grunfeld) 98pcdtest(pGrunfeld$value) 99 100## local test 101## define neighbours for individual 2: 1, 3, 4, 5 in lower triangular matrix 102w <- matrix(0, ncol= 10, nrow=10) 103w[2,1] <- w[3,2] <- w[4,2] <- w[5,2] <- 1 104pcdtest(g, w = w) 105 106 107## cortab 108pGrunfeld <- pdata.frame(Grunfeld) 109grp <- c(rep(1, 100), rep(2, 50), rep(3, 50)) # make 3 groups 110cortab(pGrunfeld$value, grouping = grp, groupnames = c("A", "B", "C")) 111 112## ercomp 113data("Produc", package = "plm") 114ercomp(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, 115 method = "walhus", effect = "time") 116z <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 117 data = Produc, random.method = "walhus", 118 effect = "time", model = "random") 119ercomp(z) 120ercomp(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, 121 method = "amemiya", effect = "twoways") 122 123## index 124data("Grunfeld", package = "plm") 125Gr <- pdata.frame(Grunfeld, index = c("firm", "year")) 126m <- plm(inv ~ value + capital, data = Gr) 127index(Gr, "firm") 128index(Gr, "time") 129index(Gr$inv, c(2, 1)) 130index(m, "id") 131 132# with additional group index 133data("Produc", package = "plm") 134pProduc <- pdata.frame(Produc, index = c("state", "year", "region")) 135index(pProduc, 3) 136index(pProduc, "region") 137index(pProduc, "group") 138 139## is.pbalanced 140Grunfeld_missing_period <- Grunfeld[-2, ] 141is.pbalanced(Grunfeld_missing_period) # check if balanced: FALSE 142pdim(Grunfeld_missing_period)$balanced # same 143pGrunfeld_missing_period <- pdata.frame(Grunfeld_missing_period) 144is.pbalanced(Grunfeld_missing_period) 145is.pbalanced(pGrunfeld_missing_period$inv) 146 147## is.pconsecutive 148is.pconsecutive(Grunfeld) 149is.pconsecutive(Grunfeld, index=c("firm", "year")) 150 151# delete 2nd row (2nd time period for first individual) 152# -> non consecutive 153Grunfeld_missing_period <- Grunfeld[-2, ] 154is.pconsecutive(Grunfeld_missing_period) 155all(is.pconsecutive(Grunfeld_missing_period)) # FALSE 156 157# delete rows 1 and 2 (1st and 2nd time period for first individual) 158# -> consecutive 159Grunfeld_missing_period_other <- Grunfeld[-c(1,2), ] 160is.pconsecutive(Grunfeld_missing_period_other) # all TRUE 161 162# delete year 1937 (3rd period) for _all_ individuals 163Grunfeld_wo_1937 <- Grunfeld[Grunfeld$year != 1937, ] 164is.pconsecutive(Grunfeld_wo_1937) # all FALSE 165 166# pdata.frame interface 167pGrunfeld <- pdata.frame(Grunfeld) 168pGrunfeld_missing_period <- pdata.frame(Grunfeld_missing_period) 169is.pconsecutive(pGrunfeld) # all TRUE 170is.pconsecutive(pGrunfeld_missing_period) # first FALSE, others TRUE 171 172# panelmodel interface (first, estimate some models) 173mod_pGrunfeld <- plm(inv ~ value + capital, data = Grunfeld) 174mod_pGrunfeld_missing_period <- plm(inv ~ value + capital, data = Grunfeld_missing_period) 175 176is.pconsecutive(mod_pGrunfeld) 177is.pconsecutive(mod_pGrunfeld_missing_period) 178 179nobs(mod_pGrunfeld) # 200 180nobs(mod_pGrunfeld_missing_period) # 199 181 182# pseries interface 183pinv <- pGrunfeld$inv 184pinv_missing_period <- pGrunfeld_missing_period$inv 185 186is.pconsecutive(pinv) 187is.pconsecutive(pinv_missing_period) 188 189# default method for arbitrary vectors or NULL 190inv <- Grunfeld$inv 191inv_missing_period <- Grunfeld_missing_period$inv 192is.pconsecutive(inv, id = Grunfeld$firm, time = Grunfeld$year) 193is.pconsecutive(inv_missing_period, id = Grunfeld_missing_period$firm, 194 time = Grunfeld_missing_period$year) 195 196# only id and time are needed for evaluation 197is.pconsecutive(NULL, id = Grunfeld$firm, time = Grunfeld$year) 198 199## is.pseries 200Em <- pdata.frame(EmplUK) 201z <- Em$output 202 203class(z) # pseries as indicated by class 204is.pseries(z) # and confirmed by check 205 206# destroy index of pseries and re-check 207attr(z, "index") <- NA 208is.pseries(z) # now FALSE 209 210## model.frame, model.matrix 211pGrunfeld <- pdata.frame(Grunfeld) 212 213# then make a model frame from a formula and a pdata.frame 214form <- inv ~ value 215mf <- model.frame(pGrunfeld, form) 216 217# then construct the (transformed) model matrix (design matrix) 218# from model frame 219modmat <- model.matrix(mf, model = "within") 220 221## retrieve model frame and model matrix from an estimated plm object 222fe_model <- plm(form, data = pGrunfeld, model = "within") 223model.frame(fe_model) 224model.matrix(fe_model) 225 226# same as constructed before 227all.equal(mf, model.frame(fe_model), check.attributes = FALSE) # TRUE 228all.equal(modmat, model.matrix(fe_model), check.attributes = FALSE) # TRUE 229 230 231## pmodel.response 232 233form <- inv ~ value + capital 234mf <- model.frame(pGrunfeld, form) 235# construct (transformed) response of the within model 236resp <- pmodel.response(form, data = mf, model = "within", effect = "individual") 237# retrieve (transformed) response directly from model frame 238resp_mf <- pmodel.response(mf, model = "within", effect = "individual") 239 240# retrieve (transformed) response from a plm object, i.e., an estimated model 241fe_model <- plm(form, data = pGrunfeld, model = "within") 242pmodel.response(fe_model) 243 244# same as constructed before 245all.equal(resp, pmodel.response(fe_model), check.attributes = FALSE) # TRUE 246 247 248 249## nobs 250z <- plm(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp,data=Produc, 251 model="random", subset = gsp > 5000) 252 253nobs(z) # total observations used in estimation 254pdim(z)$nT$N # same information 255pdim(z) # more information about the dimensions (no. of individuals and time periods) 256 257# illustrate difference between nobs and pdim for first-difference model 258data("Grunfeld", package = "plm") 259fdmod <- plm(inv ~ value + capital, data = Grunfeld, model = "fd") 260nobs(fdmod) # 190 261pdim(fdmod)$nT$N # 200 262 263## pgmm 264## Arellano and Bond (1991), table 4 col. b 265z1 <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) 266 + log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99), 267 data = EmplUK, effect = "twoways", model = "twosteps") 268summary(z1, robust = FALSE) 269 270## Blundell and Bond (1998) table 4 (cf. DPD for OX p. 12 col. 4) 271z2 <- pgmm(log(emp) ~ lag(log(emp), 1)+ lag(log(wage), 0:1) + 272 lag(log(capital), 0:1) | lag(log(emp), 2:99) + 273 lag(log(wage), 2:99) + lag(log(capital), 2:99), 274 data = EmplUK, effect = "twoways", model = "onestep", 275 transformation = "ld") 276summary(z2, robust = TRUE) 277 278# Same with the old formula or dynformula interface 279# Arellano and Bond (1991), table 4, col. b 280z1 <- pgmm(log(emp) ~ log(wage) + log(capital) + log(output), 281 lag.form = list(2,1,0,1), data = EmplUK, 282 effect = "twoways", model = "twosteps", 283 gmm.inst = ~log(emp), lag.gmm = list(c(2,99))) 284summary(z1, robust = FALSE) 285 286## Blundell and Bond (1998) table 4 (cf DPD for OX p. 12 col. 4) 287z2 <- pgmm(dynformula(log(emp) ~ log(wage) + log(capital), list(1,1,1)), 288 data = EmplUK, effect = "twoways", model = "onestep", 289 gmm.inst = ~log(emp) + log(wage) + log(capital), 290 lag.gmm = c(2,99), transformation = "ld") 291summary(z2, robust = TRUE) 292 293## pht (deprecated) 294# deprecated way with pht() for HT 295data("Wages", package = "plm") 296ht <- pht(lwage ~ wks + south + smsa + married + exp + I(exp^2) + 297 bluecol + ind + union + sex + black + ed | 298 sex + black + bluecol + south + smsa + ind, 299 data = Wages, model = "ht", index = 595) 300summary(ht) 301 302am <- pht(lwage ~ wks + south + smsa + married + exp + I(exp^2) + 303 bluecol + ind + union + sex + black + ed | 304 sex + black + bluecol + south + smsa + ind, 305 data = Wages, model = "am", index = 595) 306summary(am) 307 308## pldv 309pder.avail <- if (!requireNamespace("pder", quietly = TRUE)) FALSE else TRUE 310if(pder.avail) { 311data("Donors", package = "pder") 312pDonors <- pdata.frame(Donors, index = "id") 313modA <- pldv(donation ~ treatment + prcontr, data = pDonors, 314 model = "random", method = "bfgs") 315summary(modA) 316modB <- pldv(donation ~ treatment * prcontr - prcontr, data = pDonors, 317 model = "random", method = "bfgs") 318summary(modB) 319invisible(NULL) 320} 321 322## pwartest 323pwartest(log(emp) ~ log(wage) + log(capital), data = EmplUK) 324pwartest(log(emp) ~ log(wage) + log(capital), data = EmplUK, type = "HC3") 325 326## pwfdtest 327pwfdtest(log(emp) ~ log(wage) + log(capital), data = EmplUK) 328pwfdtest(log(emp) ~ log(wage) + log(capital), data = EmplUK, h0 = "fe") 329pwfdtest(log(emp) ~ log(wage) + log(capital), data = EmplUK, type = "HC3", h0 = "fe") 330 331mod <- plm(log(emp) ~ log(wage) + log(capital), data = EmplUK, model = "fd") 332pwfdtest(mod) 333pwfdtest(mod, h0 = "fe") 334pwfdtest(mod, type = "HC3", h0 = "fe") 335 336# pwtest 337pwtest(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc) 338pwtest(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, effect = "time") 339 340## panelmodel interface 341# first, estimate a pooling model, than compute test statistics 342form <- formula(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp) 343pool_prodc <- plm(form, data = Produc, model = "pooling") 344pwtest(pool_prodc) # == effect="individual" 345pwtest(pool_prodc, effect="time") 346