如何使用 plm 计算 R 中 gmm 模型的 BIC 和 AIC?
How to calculate BIC and AIC for a gmm model in R using plm?
我正在使用 plm
库估算 GMM 模型。我有不同的矩条件。
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')"
关于如何计算 BIC 和 AIC 或在不同矩限制之间进行选择的替代方法有什么想法吗?
我正在关注这个 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)
}
需要考虑不同版本的面板模型的不同尺寸(和参数数量)。
继续前面的示例:
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)
}
}
}
我正在使用 plm
库估算 GMM 模型。我有不同的矩条件。
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')"
关于如何计算 BIC 和 AIC 或在不同矩限制之间进行选择的替代方法有什么想法吗?
我正在关注这个 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)
}
需要考虑不同版本的面板模型的不同尺寸(和参数数量)。 继续前面的示例:
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)
}
}
}