3

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?

Helix123
  • 3,502
  • 2
  • 16
  • 36
Mario GS
  • 859
  • 8
  • 22

2 Answers2

5

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)
}
dstrants
  • 7,423
  • 2
  • 19
  • 27
Alex
  • 66
  • 1
  • 2
  • 1
    Maybe one thing is worthwile to point out? I just got stuck due to this: The effect = 'twoway' plm model and the formula by Alex will not include the time and individual effects in 'p' (the number of parameters) here. In your original question, you could write a dummy regression and then AIC() would include these dummies in 'p'. Usually you probably don't want this, though, but its still important to make sure what we compare. – Jakob Sep 03 '19 at 13:47
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)
    }
    }
    }
Rookie
  • 21
  • 2