Driscoll and Kraay standard errors in FE regression: reproducing Stata xtscc output in R

Question

I am trying to replicate the results provided by the Stata command xtscc in R with package plm but I am having some trouble to see the same standard errors I am using a dataset from the package plm also in Stata for replication purposes.

# code to obtain dataset
library(lmtest)
library(car)
library(tidyverse)
data("Produc", package="plm")
write.dta(Produc,"test.dta")

My aim is to run a two way-fixed effect panel model estimation with Driscoll and Kraay standard errors. The routine in Stata is the following

use "test.dta", clear \\ to import data
** i declare the panel 
xtset state year
* create the dummies for the time fixed effects
quietly tab year, gen(yeardum)
* run a two way fixed effect regression model with Driscoll and Kraay standard errors
xi: xtscc gsp pcap emp unemp yeardum*,fe 
* results are the following
                    Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
        pcap |  -.1769881    .265713    -0.67   0.515    -.7402745    .3862983
         emp |   40.61522   2.238392    18.14   0.000     35.87004     45.3604
       unemp |   23.59849   85.10647     0.28   0.785    -156.8192    204.0161

In R I use the following routine:

# I declare the panel
Produc <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)
# run a two way fixed effect model
femodel <- plm(gsp~pcap+emp+unemp, data=Produc,effect = "twoway", 
               index = c("iso3c","year"), model="within")
# compute Driscoll and Kraay standard errors using vcovSCC
coeftest(femodel, vcovSCC(femodel))

pcap  -0.17699    0.25476 -0.6947   0.4874    
emp   40.61522    2.14610 18.9252   <2e-16 ***
unemp 23.59849   81.59730  0.2892   0.7725

While point estimates are the same that in Stata, standard errors are different.

To check whether I am using the "wrong" small sample adjustment for standard errors, I also tryed running the coeftest with all available adjustments, but none yields the same values as xtscc.

library(purrr)
results <- map(c("HC0", "sss", "HC1", "HC2", "HC3", "HC4"),~coeftest(femodel, vcovSCC(femodel,type = .x)))
walk(results,print)
# none of the estimated standard errors is the same as xtscc

Does anyone know how I can replicate the results of Stata in R?

Likely as Stata includes an intercept in their within model and plm does not, so small sample adjustments are slightly different (finite case; irrel. in large samples). There are discussions about the differences Stata vs. plm in the small sample adjustment here or on Stackoverflow, maybe that helps. You can also look into the coding (or the documentation) of plm::within_intercept to see how the FE model with an intercept is calculated (documentation was references). I suspect xtssc somehow uses Stata's xtreg. I don't know if xtssc uses the same small sample adjustment as xtreg. — Helix123
Thanks for this comment, the problem with the real data, is that while R coefficients are significant, the Stata coefficients are not. the adjustment is of 0.03. Also is there any chanche that you know what type of standard errors are used in Stata? From the documentation of xtscc is not clear.. — Alex
xtscc is afaik a user contributed command. If its documentation does not reveal the small sample adjustment employed, you might want to read the functions source code or contact its author. In plm, type = sss mimics Stata's small sample adjustment (hence its name). — Helix123
Note that the index argument in your call to plm is ignored as your data is already a pdata.frame. Also, the index argument contains a variable not in the data (iso3c). — Helix123

Helix123 Helix123 · Accepted Answer · 2020-12-12T21:54:32

Since plm version 2.4, its function within_intercept(., return.model = TRUE) can return the full model of a within model with the intercept as in Stata. With this, it is possible to exactly replicate the result of Stata's user contributed command xtscc.

The way xtscc seems to work is by estimating the twoway FE model as a one-way FE model + dummies for the time dimension. So let's replicate that with plm:

data("Produc", package="plm")
Produc <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)

femodel <- plm(gsp ~ pcap + emp + unemp + factor(year), data = Produc, model="within")
femodelint  <- within_intercept(femodel,  return.model = TRUE)

lmtest::coeftest(femodelint, vcov. = function(x) vcovSCC(x, type = "sss"))
#                     Estimate  Std. Error t value              Pr(>|t|)    
# (Intercept)      -6547.68816  3427.47163 -1.9104             0.0564466 .  
# pcap                -0.17699     0.26571 -0.6661             0.5055481    
# emp                 40.61522     2.23839 18.1448 < 0.00000000000000022 ***
# unemp               23.59849    85.10647  0.2773             0.7816356    
# [...]

Driscoll and Kraay standard errors in FE regression: reproducing Stata xtscc output in R

1 Answers