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:

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, 
                     sd = sd_eps * sqrt(fracFirmEps)), 
               each=nrT) + 
                 rnorm(nrN * nrT, mean=0, sd = sd_eps * sqrt(1-fracFirmEps))
    X   <- rep(rnorm(nrN, 
                     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)
     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.

  • $\begingroup$ The manual states "If the argument vcov was specified, the standard errors (stored as attribute "se" in the return value) are the respective robust standard errors". I don't see that you used the vcov argument above: user11562's shows you one way to use the vcov argument. $\endgroup$
    – PatrickT
    Commented May 17, 2022 at 2:07

1 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?

  • $\begingroup$ Sorry for the late delay; anyways, your approach doesn't really change my standard errors and your proposed solution is going again in the direction of the link I provided above. That is, computing the standard errors yourself. I am just wondering why plm is not doing that automatically? $\endgroup$ Commented May 30, 2012 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.