context("plm objects - first differences model")
set.seed(20190513)
skip_if_not_installed("plm")
skip_if_not_installed("AER")
library(plm, quietly=TRUE)
data(Fatalities, package = "AER")
Fatalities <- within(Fatalities, {
frate <- 10000 * fatal / pop
drinkagec <- cut(drinkage, breaks = 18:22, include.lowest = TRUE, right = FALSE)
drinkagec <- relevel(drinkagec, ref = 4)
})
plm_FD <- plm(frate ~ beertax + drinkagec + miles + unemp + log(income),
data = Fatalities, index = c("state", "year"),
model = "fd")
n_obs <- nobs(plm_FD)
target <- with(Fatalities, 1 / pop[year != levels(year)[1]])
test_that("bread works", {
y <- na.omit(diff(plm_FD$model$frate))
cluster <- findCluster.plm(plm_FD)
expect_true(check_bread(plm_FD, cluster = cluster, y = y))
sigma_sq <- with(plm_FD, sum(residuals^2) / df.residual)
expect_equal(vcov(plm_FD), bread(plm_FD) * sigma_sq / v_scale(plm_FD))
})
test_that("CR0 and CR1S agree with arellano vcov", {
expect_equal(vcovHC(plm_FD, method="arellano", type = "HC0", cluster = "group"),
as.matrix(vcovCR(plm_FD, type = "CR0")), check.attributes = FALSE)
expect_equal(vcovHC(plm_FD, method="arellano", type = "sss", cluster = "group"),
as.matrix(vcovCR(plm_FD, type = "CR1S")), check.attributes = FALSE)
})
test_that("CR0 and CR1S agree with arellano vcov for versions <= 2.6-1", {
skip_if(packageVersion("plm") > "2.6-1")
X <- model_matrix(plm_FD)
e <- residuals(plm_FD)
index <- attr(model.frame(plm_FD), "index")
cluster <- index[[2]]
cluster <- cluster[index[[2]] != levels(index[[2]])[1]]
estmats <- sapply(split.data.frame(e * X, cluster, drop = TRUE), colSums)
meat <- tcrossprod(estmats)
bread <- chol2inv(chol(crossprod(X)))
vcov_time <- bread %*% meat %*% bread
attr(vcov_time, "dimnames") <- attr(meat, "dimnames")
expect_equal(vcov_time, as.matrix(vcovCR(plm_FD, cluster = "time", type = "CR0")))
baloney <- tcrossprod(estmats[,-6])
vcov_baloney <- bread %*% baloney %*% bread
attr(vcov_baloney, "dimnames") <- attr(baloney, "dimnames")
expect_equal(vcov_baloney,
vcovHC(plm_FD, method="arellano", type = "HC0", cluster = "time"), check.attributes = FALSE)
})
test_that("CR0 and CR1S agree with arellano vcov for versions > 2.6-1", {
skip_if_not_installed("plm", minimum_version = "2.6-2")
expect_equal(vcovHC(plm_FD, method="arellano", type = "HC0", cluster = "time"),
as.matrix(vcovCR(plm_FD, type = "CR0", cluster = "time")), check.attributes = FALSE)
expect_equal(vcovHC(plm_FD, method="arellano", type = "sss", cluster = "time"),
as.matrix(vcovCR(plm_FD, type = "CR1S", cluster = "time")), check.attributes = FALSE)
})
test_that("vcovCR options work for CR2", {
CR2_iv <- vcovCR(plm_FD, type = "CR2")
expect_equal(vcovCR(plm_FD, cluster = Fatalities$state, type = "CR2"), CR2_iv)
expect_equal(vcovCR(plm_FD, type = "CR2", inverse_var = TRUE), CR2_iv)
expect_equal(vcovCR(plm_FD, type = "CR2",
target = rep(1, n_obs),
inverse_var = TRUE), CR2_iv)
CR2_not <- vcovCR(plm_FD, type = "CR2", inverse_var = FALSE)
expect_equivalent(CR2_not, CR2_iv)
expect_equal(vcovCR(plm_FD, cluster = Fatalities$state, type = "CR2", inverse_var = FALSE), CR2_not)
expect_equal(vcovCR(plm_FD, type = "CR2", target = rep(1, n_obs)), CR2_not)
expect_equal(vcovCR(plm_FD, type = "CR2", target = rep(1, n_obs), inverse_var = FALSE), CR2_not)
expect_false(identical(vcovCR(plm_FD, type = "CR2", target = target), CR2_not))
})
test_that("vcovCR options work for CR4", {
skip_on_cran()
CR4_iv <- vcovCR(plm_FD, type = "CR4")
expect_equal(vcovCR(plm_FD, cluster = Fatalities$state, type = "CR4"), CR4_iv)
expect_equal(vcovCR(plm_FD, type = "CR4", inverse_var = TRUE), CR4_iv)
expect_equal(vcovCR(plm_FD, type = "CR4", target = rep(1, n_obs), inverse_var = TRUE), CR4_iv)
CR4_not <- vcovCR(plm_FD, type = "CR4", inverse_var = FALSE)
expect_equivalent(CR4_not, CR4_iv, tolerance = 10^-3)
expect_equal(vcovCR(plm_FD, cluster = Fatalities$state, type = "CR4", inverse_var = FALSE), CR4_not)
expect_equal(vcovCR(plm_FD, type = "CR4", target = rep(1, n_obs)), CR4_not)
expect_equal(vcovCR(plm_FD, type = "CR4", target = rep(1, n_obs), inverse_var = FALSE), CR4_not)
expect_false(identical(vcovCR(plm_FD, type = "CR4", target = target), CR4_not))
})
test_that("CR2 is target-unbiased", {
expect_true(check_CR(plm_FD, vcov = "CR2", tol = 10^-7))
expect_true(check_CR(plm_FD, vcov = "CR2", inverse_var = FALSE))
expect_true(check_CR(plm_FD, cluster = "time", vcov = "CR2", tol = 10^-7))
expect_true(check_CR(plm_FD, cluster = "time", vcov = "CR2", inverse_var = FALSE))
})
test_that("vcovCR is equivalent to vcovHC when clusters are all of size 1", {
CR_types <- paste0("CR",c(0,2))
HC_types <- paste0("HC",c(0,2))
CR_individual <- lapply(CR_types, function(t)
as.matrix(vcovCR(plm_FD, cluster = 1:nrow(Fatalities), type = t)))
HC_individual <- lapply(HC_types, function(t)
vcovHC(plm_FD, method = "white1", type = t))
expect_equal(CR_individual, HC_individual, check.attributes = FALSE)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.