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

Question

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?

Answer 1

I am following the answer to this question.

For further reference on the AIC criteria, you can look at Wikipedia.

Here is the code that should work. However, you didn't provide any reproducible model estimation. Hence, this is without validation for your case.

# Function: Calculates AIC based on an lm or plm object

AIC_adj <- function(mod){
  # Number of observations
  n.N   <- nrow(mod$model)
  # Residuals vector
  u.hat <- residuals(mod)
  # Variance estimation
  s.sq  <- log( (sum(u.hat^2)/(n.N)))
  # Number of parameters (incl. constant) + one additional for variance estimation
  p     <-  length(coef(mod)) + 1

  # Note: minus sign cancels in log likelihood
  aic <- 2*p  +  n.N * (  log(2*pi) + s.sq  + 1 ) 

  return(aic)
}

Answer 2

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)
    }
    }
    }