代码之家 › 专栏 › 技术社区 › Mario GS

如何使用plm计算R中gmm模型的BIC和AIC?

gmm plm benchmarking r

3

Mario GS · 技术社区 · 7 年前

我正在使用 plm 图书馆我有不同的力矩条件。

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)

Z ,因此

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
}

我想使用 BIC AIC ,例如

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

2 回复 | 直到 7 年前

1

5

dstrants 7 年前

我正在寻找答案 question .

有关AIC标准的更多参考,请参阅维基百科。

以下是应该有效的代码。但是,您没有提供任何可复制的模型估计。因此,这对您的案例没有验证。

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

2

Rookie 6 年前

继续上一个示例:

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