2012-05-21 17 views
6

il mio problema è questo: ottengo NA dove dovrei ottenere alcuni valori nel calcolo degli errori standard robusti.Regressione dati pannello: errori standard robusti

Sto cercando di eseguire una regressione del pannello di effetti fissi con errori standard robusti del cluster. Per questo, seguo lo Arai (2011) che a p. 3 segue Stock/ Watson (2006) (successivamente pubblicato in Econometrica, per coloro che hanno accesso). Vorrei correggere i gradi di libertà con (M/(M-1)*(N-1)/(N-K) contro il bias al ribasso poiché il mio numero di cluster è finito e ho dati sbilanciati.

Problemi simili sono stati registrati prima [1, 2] su StackOverflow e problemi correlati [3] su CrossValidated.

Arai (e la risposta nel 1 ° link) utilizza il seguente codice per le funzioni (che fornisco i miei dati qui di seguito con qualche ulteriore commento):

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

, dove il gcenter calcola le deviazioni dalla media (effetto fisso). Quindi continuo e faccio la regressione con DS_CODE come variabile del mio cluster (ho chiamato i miei dati "dati").

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

e desidera calcolare

clx(datalm, dfcw, data$DS_CODE) 

Tuttavia, quando ho voglia di calcolare uj (vedi sopra formula clx) per la varianza, ottengo solo all'inizio alcuni valori per i miei regressori, poi molti zeri. Se questo input è uj utilizzato per la varianza, solo il risultato è NAs.

I miei dati

Dato che i miei dati possono essere di struttura speciale e io non riesco a capire il problema, ho posto la cosa intera come link da Hotmail. Il motivo è che con altri dati (presi da Arai (2011)) il mio problema non si verifica. Scusa in anticipo per il casino, ma sarei molto grato se potessi comunque dare un'occhiata a questo. Il file è un file .txt 5mb contenente dati puramente.

+0

Arai non esistono più sotto il vostro link. Puoi fornire il link attuale? – MERose

risposta

2

Dopo un po 'di tempo a giocare in giro, funziona per me e mi dà:

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

Questo ci lascia con la domanda perché non lo fa per voi. Immagino che abbia qualcosa a che fare con il formato dei tuoi dati. Tutto è numerico? Ho convertito le classi della colonna e sembra che per me:

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

Che cosa significa str(data) ritorno per voi?

+0

Grazie mille per il tuo impegno e la tua risposta! Il mio 'str (data)' restituisce * Fattore * per 'DS_CODE' e * int * per' DNEW'. Tutti gli altri risultati sono uguali .... MA: Questa è la cosa più strana: funziona ora se uso il set di dati * ridotto * (ti ho dato solo il set di dati piccolo senza i miei altri varialbes e i numeri di riga R). Con il grande set, ottengo 1 singola riga di 'NAs 'nel calcolo di * uj *. Se esporto l'intero set di dati SENZA numeri di riga ('row.names = FALSE'), importalo di nuovo e esegui la regressione, funziona con il big data set. Non so perché ... – Jan

+0

Felice che funzioni ora. –

0

Il pacchetto plm può stimare SE in cluster per regressioni di pannello. I dati originali non sono più disponibili, quindi ecco un esempio utilizzando i dati fittizi.

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 

Se stai usando lm modelli (invece di plm), poi il pacchetto multiwayvcov può aiutare.

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 

Per maggiori dettagli si veda:

Consulta anche: carta

Problemi correlati