2012-05-21 55 views
6

我的问题是这样的:我得到NA我应该在计算可靠的标准错误时得到一些值。面板数据回归:稳健的标准错误

我正在尝试使用群集健全标准错误进行固定效果面板回归。为此,我按照Arai (2011)谁在页。 3遵循Stock/ Watson (2006)(稍后发布在Econometrica,对于那些有访问权限的人)。我想通过(M/(M-1)*(N-1)/(N-K)来纠正自由度,因为我的群集数量是有限的,而且我有不平衡的数据。

类似问题[12]在计算器上和交叉验证相关问题[3]之前已经公布。

新井(以及在第一环节的答案)使用的功能,下面的代码(我在下面提供了一些进一步的评论我的数据):

gcenter <- function(df1,group) { 
    variables <- paste(
     rep("C", ncol(df1)), colnames(df1), sep=".") 
    copydf <- df1 
    for (i in 1:ncol(df1)) { 
     copydf[,i] <- df1[,i] - ave(df1[,i], group,FUN=mean)} 
    colnames(copydf) <- variables 
    return(cbind(df1,copydf))} 

# 1-way adjusting for clusters 
clx <- function(fm, dfcw, cluster){ 
    # R-codes (www.r-project.org) for computing 
    # clustered-standard errors. Mahmood Arai, Jan 26, 2008. 
    # The arguments of the function are: 
    # fitted model, cluster1 and cluster2 
    # You need to install libraries `sandwich' and `lmtest' 
    # reweighting the var-cov matrix for the within model 
    library(sandwich);library(lmtest) 
    M <- length(unique(cluster)) 
    N <- length(cluster)   
    K <- fm$rank       
    dfc <- (M/(M-1))*((N-1)/(N-K)) 
    uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum)); 
    vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw 
    coeftest(fm, vcovCL) } 

,其中gcenter计算偏离平均值(固定效果)。然后,我继续并使用DS_CODE作为我的群集变量(我已将其数据命名为“数据”)进行回归。

centerdata <- gcenter(data, data$DS_CODE) 
datalm <- lm(C.L1.retE1M ~ C.MCAP_SEC + C.Impact_change + C.Mom + C.BM + C.PD + C.CashGen + C.NITA + C.PE + C.PEdummy + factor(DS_CODE), data=centerdata) 
M <- length(unique(data$DS_CODE)) 
dfcw <- datalm$df/(datalm$df - (M-1)) 

,并要计算

clx(datalm, dfcw, data$DS_CODE) 

然而,当我想计算UJ(见上述公式clx)的变化,我只得到一开始对我的回归系数的一些值,然后很多零。如果此输入uj用于方差,则只有NAs结果。

我的数据

由于我的数据可能是特殊的结构,我无法弄清楚这个问题,我发布了整个事情从Hotmail一个link。原因是有了其他数据(摘自Arai(2011)),我的问题不会发生。对不起,我很抱歉,如果你可以看看它,但我会很感激。 该文件是一个纯粹包含数据的5mb .txt文件。

+0

新居的论文你的链接下不复存在。你能提供实际的链接吗? – MERose

回答

2

一段时间玩耍后,它为我工作,给我:

      Estimate Std. Error t value Pr(>|t|)  
(Intercept)   4.5099e-16 5.2381e-16 0.8610 0.389254  
C.MCAP_SEC   -5.9769e-07 1.2677e-07 -4.7149 2.425e-06 *** 
C.Impact_change  -5.3908e-04 7.5601e-05 -7.1306 1.014e-12 *** 
C.Mom     3.7560e-04 3.3378e-03 0.1125 0.910406  
C.BM     -1.6438e-04 1.7368e-05 -9.4645 < 2.2e-16 *** 
C.PD     6.2153e-02 3.8766e-02 1.6033 0.108885  
C.CashGen    -2.7876e-04 1.4031e-02 -0.0199 0.984149  
C.NITA    -8.1792e-02 3.2153e-02 -2.5438 0.010969 * 
C.PE     -6.6170e-06 4.0138e-06 -1.6485 0.099248 . 
C.PEdummy    1.3143e-02 4.8864e-03 2.6897 0.007154 ** 
factor(DS_CODE)130324 -5.2497e-16 5.2683e-16 -0.9965 0.319028  
factor(DS_CODE)130409 -4.0276e-16 5.2384e-16 -0.7689 0.441986  
factor(DS_CODE)130775 -4.4113e-16 5.2424e-16 -0.8415 0.400089 
... 

这给我们留下了为什么它不适合你的问题。我想这与你的数据格式有关。一切都是数字吗?我转换的列类,它看起来像我:

str(dat) 
'data.frame': 48251 obs. of 12 variables: 
$ DS_CODE  : chr "902172" "902172" "902172" "902172" ... 
$ DNEW   : num 2e+05 2e+05 2e+05 2e+05 2e+05 ... 
$ MCAP_SEC  : num 78122 71421 81907 80010 82462 ... 
$ NITA   : num 0.135 0.135 0.135 0.135 0.135 ... 
$ CashGen  : num 0.198 0.198 0.198 0.198 0.198 ... 
$ BM   : num 0.1074 0.1108 0.097 0.0968 0.0899 ... 
$ PE   : num 57 55.3 63.1 63.2 68 ... 
$ PEdummy  : num 0 0 0 0 0 0 0 0 0 0 ... 
$ L1.retE1M : num -0.72492 0.13177 0.00122 0.07214 -0.07332 ... 
$ Mom   : num 0 0 0 0 0 ... 
$ PD   : num 5.41e-54 1.51e-66 3.16e-80 2.87e-79 4.39e-89 ... 
$ Impact_change: num 0 -10.59 -10.43 0.7 -6.97 ... 

什么str(data)回报吗?

+0

非常感谢您的努力和您的回答!我的'str(data)'返回* DS_CODE的* Factor *和'DNEW'的* int *。所有其他结果是相同的....但:这是最奇怪的事情:现在,如果我使用*减少*数据集(我给你只有小数据集没有我的其他varials和R行号)现在的作品。有了这个大集合,我在* uj *的计算中得到了1个单独的'NAs'。如果我没有行号('row.names = FALSE')导出我的整个数据集,再次导入并进行回归,它可以处理大数据集。我不知道为什么... – Jan

+0

现在很高兴它的作品。 –

0

plm包可以估计面板回归的聚集SE。原始数据不再可用,所以这里是一个使用虚拟数据的例子。

require(foreign) 
require(plm) 
require(lmtest) 
test <- read.dta("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta") 

fpm <- plm(y ~ x, test, model='pooling', index=c('firmid', 'year')) 

##Arellano clustered by *group* SEs 
> coeftest(fpm, vcov=function(x) vcovHC(x, cluster="group", type="HC0")) 

t test of coefficients: 

      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 0.029680 0.066939 0.4434 0.6575  
x   1.034833 0.050540 20.4755 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

如果您使用lm模型(而不是plm),那么multiwayvcov包可能会有帮助。

library("lmtest") 
library("multiwayvcov") 

data(petersen) 
m1 <- lm(y ~ x, data = petersen) 

> coeftest(m1, vcov=function(x) cluster.vcov(x, petersen[ , c("firmid")], 
    df_correction=FALSE)) 

t test of coefficients: 

      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 0.029680 0.066939 0.4434 0.6575  
x   1.034833 0.050540 20.4755 <2e-16 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

有关详细信息,请参阅:

参见: