## Replication code for: ## Robust Standard Error Estimators for Panel Models: A Unifying Approach ## 4.3. Unbalanced panels a <- array(1, dim = c(3, 3, 3)) a[1, 1, 1] <- NA apply(a, 1:2, mean) apply(a, 1:2, mean, na.rm = TRUE) ## 5.1. plm and generic sandwich estimators library("plm") data("Produc", package = "plm") fm <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp lmmod <- lm(fm, Produc) library("lmtest") library("sandwich") coeftest(lmmod, vcov = vcovHC) plmmod <- plm(fm, Produc, model = "pooling") coeftest(plmmod, vcov = vcovHC) coeftest(plmmod, vcov = function(x) vcovHC(x, method = "white1", type = "HC3")) ## 5.2. The new modular framework vcovG(plmmod, cluster = "group", inner = "cluster", l = 0) vcovHC(plmmod) vcovSCC(plmmod) myvcovDC <- function(x, ...) { Vcx <- vcovHC(x, cluster = "group", method = "arellano", ...) Vct <- vcovHC(x, cluster = "time", method = "arellano", ...) Vw <- vcovHC(x, method = "white1", ...) return(Vcx + Vct - Vw) } myvcovDC(plmmod) vcovDC(plmmod) myvcovDCS <- function(x, maxlag = NULL, ...) { w1 <- function(j, maxlag) 1 VsccL.1 <- vcovSCC(x, maxlag = maxlag, wj = w1, ...) Vcx <- vcovHC(x, cluster = "group", method = "arellano", ...) VnwL.1 <- vcovSCC(x, maxlag = maxlag, inner = "white", wj = w1, ...) return(VsccL.1 + Vcx - VnwL.1) } myvcovDCS(plmmod, maxlag = 4) ## 6. Applied examples Vw <- function(x) vcovHC(x, method = "white1") Vcx <- function(x) vcovHC(x, cluster = "group", method = "arellano") Vct <- function(x) vcovHC(x, cluster = "time", method = "arellano") Vcxt <- function(x) Vcx(x) + Vct(x) - Vw(x) Vct.L <- function(x) vcovSCC(x, wj = function(j, maxlag) 1) Vnw.L <- function(x) vcovNW(x) Vscc.L <- function(x) vcovSCC(x) Vcxt.L <- function(x) Vct.L(x) + Vcx(x) - vcovNW(x, wj = function(j, maxlag) 1) vcovs <- c(vcov, Vw, Vcx, Vct, Vcxt, Vct.L, Vnw.L, Vscc.L, Vcxt.L) names(vcovs) <- c("OLS", "Vw", "Vcx", "Vct", "Vcxt", "Vct.L", "Vnw.L", "Vscc.L", "Vcxt.L") cfrtab <- matrix(nrow = length(coef(plmmod)), ncol = 1 + length(vcovs)) dimnames(cfrtab) <- list(names(coef(plmmod)), c("Coefficient", paste("s.e.", names(vcovs)))) cfrtab[,1] <- coef(plmmod) for (i in 1:length(vcovs)) { cfrtab[, 1 + i] <- coeftest(plmmod, vcov = vcovs[[i]])[, 2] } print(t(round(cfrtab, 4))) ## 6.1. PPP regression data("Parity", package = "plm") fm <- ls ~ ld pppmod <- plm(fm, data = Parity, effect = "twoways") library("car") linearHypothesis(pppmod, "ld = 1", vcov = Vcxt.L) ## 6.2. Petersen’s artificial data petersen <- read.table(file = "test_data.txt") colnames(petersen) <- c("firmid", "year", "x", "y") ptrmod <- plm(y ~ x, data = petersen, index = c("firmid", "year"), model = "pooling") vcovs <- c(vcov, Vw, Vcx, Vct, Vcxt) names(vcovs) <- c("OLS", "Vw", "Vcx", "Vct", "Vcxt") cfrtab <- matrix(nrow = length(coef(ptrmod)), ncol = 1 + length(vcovs)) dimnames(cfrtab) <- list(names(coef(ptrmod)), c("Coefficient", paste("s.e.", names(vcovs)))) cfrtab[, 1] <- coef(ptrmod) for (i in 1:length(vcovs)) { cfrtab[, 1 + i] <- coeftest(ptrmod, vcov = vcovs[[i]])[, 2] } print(t(round(cfrtab, 4))) coeftest(ptrmod, vcov = function(x) vcovDC(x, type = "sss"))[, 2]