Solved – Newey West Covariance in SAS and in R are different

neweywestrregressionsas

I am using the proc autoreg procedure in SAS and the sandwich library in R to calculate Newey West covariance for a linear model. The two methods give the same coefficient estimates (expected because they are OLS estimators), but different t-statistics hence p-values. Below is my code for model summary in R:

library(sandwich)
model = lm(response ~ var1 + var2 + var3 + var4 + var5 + var6,data=SAS_Q)
coeftest(model, vocv=NeweyWest(model,lag=2, prewhite=FALSE))

Below is my code to generate model summary in SAS:

proc autoreg data= SAS_Q plots=all;
    model response = var1 var2 var3 var4 var5 var6/covest=neweywest;;
    ods output ParameterEstimates = _auto_full_params_;
    output out=Autoreg_Out p=FORECAST_t rm=r pm=pm lcl=lcl lclm=lclm ucl=ucl uclm=uclm;
run;
quit;

No matter how I assign the lag, I can not get the R result to match that of the SAS. More specifically, the p-value for var5 is 0.11 in R but 0.03 in SAS. It also seems that my R-output is quite immune to the change of lag.

Am I doing anything wrong in R? If so, what is it?

Best Answer

vocv should be vcov - R cannot take that part of the code into account. I would have expected R to throw an error so that you do not get a result at all, or at least a warning, but that does not seem to be the case: 

library(sandwich)
library(lmtest) # note this is needed for coeftest
y <- rnorm(100)
x <- rnorm(100)
model <- lm(y ~ x)
coeftest(model) # default OLS setting
coeftest(model, vocv=NeweyWest(model,lag=2, prewhite=FALSE)) # with typo - same results as default
coeftest(model, vcov=NeweyWest(model,lag=2, prewhite=FALSE)) # without typo - different results (only slightly so here as my data does not have autocorrelation so that HAC standard erros are not necessary

returns

    > coeftest(model)

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0687659  0.1041773 -0.6601   0.5107
x            0.0068797  0.0991511  0.0694   0.9448

> coeftest(model, vocv=NeweyWest(model,lag=2, prewhite=FALSE))

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0687659  0.1041773 -0.6601   0.5107
x            0.0068797  0.0991511  0.0694   0.9448

> coeftest(model, vcov=NeweyWest(model,lag=2, prewhite=FALSE))

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0687659  0.1075861 -0.6392   0.5242
x            0.0068797  0.0808703  0.0851   0.9324

I am not sufficiently skilled in R to give an authoritative answer (nor am I sure it is on topic for CV), but the underlying issue seems to be the combined presence of default values for existing arguments combined with ... allowing for additional arguments to be passed on.

Consider

ExampleFuncWithAdditionalArguments <- function(xy=0,...){
  xy^2
}
ExampleFuncWithAdditionalArguments(xy=5)
ExampleFuncWithAdditionalArguments(xz=5)

which returns

> ExampleFuncWithAdditionalArguments(xy=5)
[1] 25
> ExampleFuncWithAdditionalArguments(xz=5)
[1] 0

The first call is as intended, the second has a typo. R then interprets the argument xz as an additional argument, with which the function however does nothing and ExampleFuncWithAdditionalArguments falls back to the default xy=0.