Replace std.error of OLS by HAC estimator, and use stargazer to print table

It is my first time posting questions here.

I'm dealing with a case where I use time series data to do simple ols estimate, by using dynlm. Also, I want to replace the standard error by HAC estimator with some conventional truncated value as the lag variable. The code is like what follows:

dat <- ts(data)
reg <- dynlm(y ~ x1 + x2 + x3, data = dat)
ols <- summary(reg)
robust <- coeftest(reg, vcov=NeweyWest(reg, lag = round(0.75 * length(time(reg))^(1/3))))
ols$coefficients[,2:4] <- robust[,2:4]
ols

I use the same method over 16 data sets, and now I want to print the results into a latex table by using stargazer. However, the package can only generate tables by the results of a linear model rather than its summary (take the preceding code as an example, stargazer can be applied to reg rather than ols), which I already know.

It seems to me that there are two ways to deal with it

Find a way to generate latex tables by using summaries of regressions.
Replace the std errors, t stats and p-values in reg rather than those in summary(reg).

It would be so nice of you to give me an aid! Thanks!

Solution

I am not sure I understand your question correctly, but if I do the answer is quite easy and has been provided in this thread: Extending Stargazer to multiwaycov

You simply have to supply the robust standard errors and p-values to the stargazer function directly with

stargazer(ols, ..., se = robust[,2], p = robust[,4], p.auto = F)

No need to supply t-statistics as they are calculated by stargazer based on the supplied std. errors and coefficients.

I would strongly suggest to read the stargazer documentation and search SO more carefully before posting a question.