Regresja danych panelu: Solidne błędy standardowe

Question

Mój problem polega na tym, że dostaję NA, gdzie powinienem otrzymać pewne wartości w obliczeniach solidnych błędów standardowych.

Próbuję wykonać regresję z panelem efektów stałych za pomocą standardowych błędów klastrowych. W tym celu podążam za Arai (2011), który na s. 3 następuje po Stock/ Watson (2006) (później opublikowane w Econometrica, dla tych, którzy mają dostęp). Chciałbym skorygować stopnie swobody przez (M/(M-1)*(N-1)/(N-K) przeciwko odchyleniu w dół, ponieważ moja liczba klastrów jest skończona i mam niezbalansowane dane.

Podobne problemy zostały wysłane przed [1, 2] na StackOverflow i związane z tym problemy [3] na CrossValidated.

Arai (a odpowiedź w 1 linku) wykorzystuje następujący kod dla funkcji (I dostarczyć swoje dane poniżej pewnego dalszego komentarza):

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

, gdzie gcenter oblicza odchylenia od średniej (poprawiony efekt). Następnie kontynuuję i wykonuję regres z DS_CODE będącą moją zmienną klastra (nazwa moich danych to "dane").

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

i chcemy obliczyć

clx(datalm, dfcw, data$DS_CODE) 

Jednak, gdy chcę, aby obliczyć UJ (patrz wzór clx powyżej) dla wariancji, mam tylko na początku pewne wartości dla moich regresorów, następnie wiele zer. Jeśli to wejście uj zostanie użyte dla wariancji, wyniknie tylko NAs.

moje dane

Ponieważ moje dane mogą być z specjalnej konstrukcji i nie mogę dowiedzieć się problem, mogę napisać całą rzecz jako link z Hotmail. Powodem jest to, że z innymi danymi (zaczerpniętymi z Arai (2011)) mój problem nie występuje. Przykro mi z góry za bałagan, ale byłabym bardzo wdzięczna, gdybyś mogła rzucić na to okiem. Plik to plik 5tb .txt zawierający wyłącznie dane.

Answer 1

Po pewnym czasie zabawy, to działa na mnie i daje mi:

      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 
... 

To pozostawia nas z pytaniem, dlaczego nie dla ciebie. Myślę, że ma to coś wspólnego z formatem twoich danych. Czy wszystko jest numeryczne? I konwertowane klasy kolumn i wygląda na to, że dla mnie:

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 ... 

Co str(data) zwrot dla ciebie?

Answer 2

Pakiet plm może oszacować zagregowane wartości SE dla regresji panelu. Oryginalne dane nie są już dostępne, więc oto przykład użycia danych fikcyjnych.

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 

---

Jeśli używasz lm modeli (zamiast plm), wówczas pakiet multiwayvcov może pomóc.

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 

Aby uzyskać więcej informacji zobacz:

• Fama-MacBeth and Cluster-Robust (by Firm and Time) Standard Errors in R.

Zobacz także: Papier

• Double clustered standard errors for panel data