title: Panel data regression: Robust standard errors
body: <p>my problem is this: I get <code>NA</code> where I should get some values in the computation of robust standard errors. </p>
  
  <p>I am trying to do a fixed effect panel regression with cluster-robust standard errors. For this, I follow <a href="http://people.su.se/~ma/clustering.pdf" rel="nofollow noreferrer">Arai (2011)</a> who on p. 3 follows <a href="https://www.princeton.edu/~mwatson/papers/panelse_5.pdf" rel="nofollow noreferrer">Stock/ Watson (2006)</a> (later published in <a href="http://www.jstor.org/discover/10.2307/4502057?uid=3737760&amp;uid=2&amp;uid=4&amp;sid=56192307833" rel="nofollow noreferrer">Econometrica</a>, for those who have access). I would like to correct the degrees of freedom by <code>(M/(M-1)*(N-1)/(N-K)</code>   against downward bias as my number of clusters is finite and I have unbalanced data. </p>
  
  <p>Similar problems have been posted before [<a href="https://stackoverflow.com/questions/9260766/robust-standard-errors-in-ggplot2/9264761#9264761">1</a>, <a href="https://stackoverflow.com/questions/8389843/clustered-standard-errors-r-panel-data">2</a>] on StackOverflow and related problems [<a href="https://stats.stackexchange.com/questions/28492/understanding-the-plm-package-why-arent-my-standard-errors-correct">3</a>] on CrossValidated. </p>
  
  <p>Arai (and the answer in the 1st link) uses the following code for functions (<em>I provide my data below with some further comment</em>):</p>
  
  <pre><code>gcenter &lt;- function(df1,group) {
      variables &lt;- paste(
          rep("C", ncol(df1)), colnames(df1), sep=".")
      copydf &lt;- df1
      for (i in 1:ncol(df1)) {
          copydf[,i] &lt;- df1[,i] - ave(df1[,i], group,FUN=mean)}
      colnames(copydf) &lt;- variables
      return(cbind(df1,copydf))}
  
  # 1-way adjusting for clusters
  clx &lt;- 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 &lt;- length(unique(cluster))   
      N &lt;- length(cluster)           
      K &lt;- fm$rank                        
      dfc &lt;- (M/(M-1))*((N-1)/(N-K))  
      uj  &lt;- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
      vcovCL &lt;- dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
      coeftest(fm, vcovCL) }
  </code></pre>
  
  <p>,where the <code>gcenter</code> computes deviations from the mean (fixed effect). I then continue and do the regression with <code>DS_CODE</code>being my cluster variable (I have named my data 'data').</p>
  
  <pre><code>centerdata &lt;- gcenter(data, data$DS_CODE)
  datalm &lt;- 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 &lt;- length(unique(data$DS_CODE))
  dfcw &lt;- datalm$df / (datalm$df - (M-1))
  </code></pre>
  
  <p>and want to calculate </p>
  
  <pre><code>clx(datalm, dfcw, data$DS_CODE)
  </code></pre>
  
  <p>However, when I want to compute <em>uj</em> (see formula <code>clx</code> above) for the variance, I get only at the beginning some values for my regressors, then lots of zeros. If this input <em>uj</em> is used for the variance, only <code>NAs</code> result.</p>
  
  <p><strong>My data</strong></p>
  
  <p>Since my data may be of special structure and I can't figure out the problem, I post the entire thing as a <a href="https://skydrive.live.com/redir?resid=58FBD7EAB94B2963!223&amp;authkey=!APTXLzboXcf67Hs" rel="nofollow noreferrer">link</a> from Hotmail. The reason is that with other data (taken from Arai (2011)) my problem does not occur. Sorry in advance for the mess but I'd be very grateful if you could have a look at it nevertheless. 
  The file is a 5mb .txt file containing purely data. </p>