time series regression by group with newey-West standard errors plus getting R2

Question

I have a data.table with annual returns and 2 explanatory variables for 25 different equity portfolios. I would like to estimate the same lm model for each of 25 portfolios where the standard errors are NeweyWest-corrected. So far, I'm runnning the model on each portfolio with group_by from dplyr and then correct the standard errors with coeftest from the lmtest package plus NeweyWest from the sandwich package, summarized with tidy from the broom package:

library(dplyr) 
library(broom)
library(lmtest)
library(sandwich)

regressions <- data %>%
      group_by(Portfolio) %>%
      do({fit = lm(Portfolio_return ~ x1 + x2, data = .)
      tidy(coeftest(fit, vcov. = NeweyWest(fit, prewhite = FALSE)))
      })

The questions I have:

The code gives me the coefficients plus p-values, but how do I get all the other summary stats like r2, F-test etc for all the models after the NeweyWest adjustment? I like tidy, glance and augment from from the broom package but when I run regressions %>% glance(fit) I get fatal error and R crashes.
How can I summarize all the stats (like coefficients, p-values, R2) for all 25 models in one data.table or data.frame?

Thanks a lot!

StupidWolf StupidWolf · Accepted Answer · 2020-10-03T17:12:36

I think what you need is possible now with dplyr, tidyr and purrr, below I used an example dataset:

library(dplyr) 
library(broom)
library(lmtest)
library(sandwich)
library(tidyr)
library(purrr)

data(Investment)

set.seed(999)
data = rbind(
       data.frame(Investment + matrix(rnorm(140),20,7),set="A"),
       data.frame(Investment + matrix(rnorm(140),20,7),set="B")
       )

We do something similar to do() in dplyr, just that we nest the data into a data column , fit and store the results, then apply glance and tidy:

regressions <-data %>% 
nest(data=c("GNP","Investment","Price","Interest","RealGNP","RealInv","RealInt")) %>% 
mutate(fit=map(data,~lm(RealInv ~ RealGNP + RealInt, data = .))) %>%
mutate(
glance = map(fit,glance),
tidy = map(fit,~tidy(coeftest(., vcov. = NeweyWest(., prewhite = FALSE))))
)

You can extract the tidy and glance results

regressions %>% unnest(tidy)
# A tibble: 6 x 9
  set   data      fit    glance     term   estimate std.error statistic  p.value
  <fct> <list>    <list> <list>     <chr>     <dbl>     <dbl>     <dbl>    <dbl>
1 A     <tibble … <lm>   <tibble [… (Inte…  -15.8     18.1       -0.871  3.97e-1
2 A     <tibble … <lm>   <tibble [… RealG…    0.172    0.0174     9.88   3.24e-8
3 A     <tibble … <lm>   <tibble [… RealI…   -2.55     3.25      -0.785  4.44e-1
4 B     <tibble … <lm>   <tibble [… (Inte…   -5.03    19.7       -0.256  8.01e-1
5 B     <tibble … <lm>   <tibble [… RealG…    0.161    0.0185     8.71   1.80e-7
6 B     <tibble … <lm>   <tibble [… RealI…    0.662    2.40       0.275  7.86e-1

regressions %>% unnest(glance)
# A tibble: 2 x 15
  set   data  fit   tidy  r.squared adj.r.squared sigma statistic p.value    df
  <fct> <lis> <lis> <lis>     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>
1 A     <tib… <lm>  <tib…     0.830         0.809  16.5      39.0 7.02e-7     3
2 B     <tib… <lm>  <tib…     0.801         0.776  17.6      32.2 2.47e-6     3
# … with 5 more variables: logLik <dbl>, AIC <dbl>, BIC <dbl>, deviance <dbl>,
#   df.residual <int>

time series regression by group with newey-West standard errors plus getting R2

1 Answers