How to calculate BIC and AIC for a gmm model in R using plm?

user_4304635 0
Viewed 3177

I'm estimating a GMM model using the plm library. I have different moment conditions.

Z <- list(~YDWPP + ST_DEGREE, ~YDWPP + ST_DEGREE, ~YDWPP + ST_DEGREE, 
    ~YDWPP + ST_DEGREE, ~YDWPP + ST_TRANSITIVITY, ~YDWPP + ST_STRUC_HOLE, 
    ~YDWPP + ST_STRUC_HOLE, ~YDWPP + ST_STRUC_HOLE, ~YDWPP + 
        ST_STRUC_HOLE)

Z <- lapply(Z, as.formula)

lg.gmm <- list(c(4L, 8L), c(5L, 8L), c(6L, 8L), 7:8, 7:8, c(4L, 8L), c(5L, 
8L), c(6L, 8L), 7:8)

I am running a loop for each set of moment restrictions Z, such that

out.1 <- list()
for(i in seq_along(Z)){
  plm.gmm <-
  pgmm(
  dynformula(as.formula(model), lg),
  data = pdata,
  effect = 'twoway',
  model = 'twostep',
  transformation = 'd',
  gmm.inst = Z[[i]],
  lag.gmm =  c(lg.gmm[[i]][[1]], lg.gmm[[i]][[2]])
  )
sum <- summary(plm.gmm, robust = T)
print(sum)
out.1[[i]] <- sum
}

I would like to compare these models using BIC and AIC, for instance

AIC(plm.gmm, k=2)
Error in UseMethod("logLik") : 
  no applicable method for 'logLik' applied to an object of class "c('pgmm', 'panelmodel')"

Any ideas on how to compute BIC and AIC or alternative methods to choose between different moment restrictions?

plm gmm benchmarking r
2 Answers

One need to take in consider different dimensions (and number of parameters) of various versions of panel models. Continuing from the previous example:

    aicbic_plm <- function(object, criterion) {


    # object is "plm", "panelmodel" 
    # Lets panel data has index :index = c("Country", "Time")

    sp = summary(object)

    if(class(object)[1]=="plm"){
    u.hat <- residuals(sp) # extract residuals
    df <- cbind(as.vector(u.hat), attr(u.hat, "index"))
    names(df) <- c("resid", "Country", "Time")
    c = length(levels(df$Country)) # extract country dimension 
    t = length(levels(df$Time)) # extract time dimension 
    np = length(sp$coefficients[,1]) # number of parameters
    n.N = nrow(sp$model) # number of data
    s.sq  <- log( (sum(u.hat^2)/(n.N))) # log sum of squares

    # effect = c("individual", "time", "twoways", "nested"),
    # model = c("within", "random", "ht", "between", "pooling", "fd")

    # I am making example only with some of the versions:

    if (sp$args$model == "within" & sp$args$effect == "individual"){
    n = c
    np = np+n+1 # update number of parameters
    }

    if (sp$args$model == "within" & sp$args$effect == "time"){
    T = t
    np = np+T+1 # update number of parameters
    }

    if (sp$args$model == "within" & sp$args$effect == "twoways"){
    n = c
    T = t
    np = np+n+T # update number of parameters
    }
    aic <- round(       2*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)
    bic <- round(log(n.N)*np  +  n.N * (  log(2*pi) + s.sq  + 1 ),1)

    if(criterion=="AIC"){
    names(aic) = "AIC"
    return(aic)
    }
    if(criterion=="BIC"){
    names(bic) = "BIC"
    return(bic)
    }
    }
    }
user_25251570
Related