Understanding the plm package - Why aren't my standard errors correct?

Question

For fun, I tried to replicate the results of Petersen (2009) who deals with the correct estimation of standard errors in finance panel data sets.

In a nutshell, he estimates the following standard regression for a panel data set:

$$ Y_{it} = X_{it} \beta + \epsilon_{it} $$

where $\epsilon_{it} = \gamma_i + \eta_{it}$ and $x_{it} = \mu_{i} + \nu_{it}$. Hence, both the residual and the independent variable have a firm-specific component. Petersen goes on to show that this results in biased standard errors when applying the standard OLS. For example, he shows in table 1 of his paper that if both the residual volatility and the variable volatility are driven by 50% by a firm-specific component, the true standard errors are nearly twice as large as the ones given by OLS.

He shows that in a MCS and I reproduced those results in R, as you can see from the code below. Naturally, I asked myself how I would compute the correct standard errors in R and the package of choice seemed to be plm. However, I just don't get the correct results out of it and I don't know what I miss.

Here is my code:

library(plm)
runMCS <- function(runs, nrN, nrT, fracFirmX, fracFirmEps, sd_X, sd_eps, beta) {

  betas    <- numeric(runs)
  se_betas <- numeric(runs)
  panel_betas    <- numeric(runs)
  se_panel_betas <- numeric(runs)

  for (i in 1:runs) {

    #Model epsilon, X, and Y
    eps <- rep(rnorm(nrN, 
                     mean=0, 
                     sd = sd_eps * sqrt(fracFirmEps)), 
               each=nrT) + 
                 rnorm(nrN * nrT, mean=0, sd = sd_eps * sqrt(1-fracFirmEps))
    X   <- rep(rnorm(nrN, 
                     mean=0, 
                     sd = sd_X   * sqrt(fracFirmX)),   
               each=nrT) + 
                 rnorm(nrN * nrT, mean=0, sd = sd_X   * sqrt(1-fracFirmX))
    Y   <- beta * X + eps

    #Compute regression (OLS)
    reg <- summary(lm(Y ~ X))

    #Save results
    betas[i]    <- reg$coef[2, 1]
    se_betas[i] <- reg$coef[2, 2]

    #Try plm
    df <- data.frame(Firm = rep(1:nrN, each=nrT),
                     Time = rep(1:nrT, times=nrN),
                     Y = Y,
                     X = X)
    preg <- summary(plm(Y ~ X, data=df, index=c("Firm", "Time"), effect="individual", model="within")) #within is fixed effects
    panel_betas[i]    <- preg$coef[1, 1]
    se_panel_betas[i] <- preg$coef[1, 2]
  }

  return(c(avg_beta = mean(betas), 
           true_se = sd(betas), 
           avg_se = mean(se_betas), 
           avg_clustered = mean(panel_betas),
           se_clustered = mean(se_panel_betas)))

}
MCS_50_50 <- runMCS(50, 500, 10, 0.5, 0.5, 1, 2, 1)
MCS_50_50
     avg_beta       true_se        avg_se avg_clustered  se_clustered 
   1.00503955    0.06020203    0.02825567    1.00433092    0.02985546

Note that I only run the simulation 50 times here because the plm function slows it down considerably. So basically, it makes virtually no difference if I call lm or plm. I'm pretty confident that I set the index and model option correct after reading the vignette of the package. However, I must miss something here! Interestingly, the package also has the fixef function and if I call that on one run, I get something like this:

summary(fixef(plm(Y ~ X, data=df, index=c("Firm", "Time"), effect="individual", model="within")))
1      13.60377     0.44112    30.8391 < 2.2e-16 ***
2    -830.74707     0.44136 -1882.2236 < 2.2e-16 ***
3    -326.96042     0.44137  -740.7840 < 2.2e-16 ***
4     169.16463     0.44246   382.3287 < 2.2e-16 ***
...

I'm not quite sure how to interpret those results, but here, I get considerably larger standard errors for each firm separately. If I would average those, I would end up with something above 0.44 which is considerably closer to the true standard errors, but still not right.

So, again a very long question from me, sorry for that ;-) Note that I did check answers before and I found this interesting link. The white paper that is referred to in the answer is interestingly the same person that implemented the solution on Petersen's webpage. So I'm pretty sure that I could get the correct standard errors by implementing Mahmood Arai's solution. But I'm looking for an already implemented and therefore safe option and I just wonder why that plm function does not work.

Answer 1

Isn't it just

coeftest(xxx, vcov=function(x) vcovHC(x, cluster=c("group","time"))) 

where xxx is the outcome of the plm() function?